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I.  EXECUTIVE  SUMMARY 

In  the  framework  of  the  present  Grant  the  scope  of  past  investigations  on  instability  of  cavity  flows  has 
been  broadened  and  unified  under  the  point  of  view  of  (nonparallel)  BiGlobal  linear  instability  theory.  The 
understanding  and  classification  of  hydrodynamic  (shear-layer-  and  in  particular  wake-mode)  and  aeroacoustic 
instabilities  has  been  advanced  by  the  first  application  of  the  author’s  residuals  algorithm  to  compressible 
open  cavity  flows.  It  has  been  shown  that  all  known  classes  of  linear  instabilities  of  open  cavity  flows,  as  well 
as  previously  unknown  linear  disturbances  in  a  cavity,  may  be  encountered  as  particular  features  of  the  same 
BiGlobal  eigenvector  at  appropriate  parameter  ranges. 

The  role  of  three-dimensionality  on  BiGlobal  cavity  flow  instability  has  been  examined;  by  reference  to 
a  theoretically  and  experimentally  well-documented  model  problem  it  has  been  shown  that  linearly  unstable 
eigendisturbances  of  the  two-dimensional  steady  basic  flow  persist  when  that  basic  state  becomes  three- 
dimensional.  This  provided  motivation  for  the  development  of  the  first  direct  solution  of  the  BiGlobal 
eigenvalue  problem  in  compressible  flow.  Accurate  and  efficient  numerical  tools  were  subsequently 
developed  in  order  to  address  three-dimensional  instability  of  both  the  classic  (empty)  as  well  as  full-bay 
models  of  cavity  flows.  New  linear  instability  mechanisms  have  been  discovered  and  independently 
confirmed  by  a  linearized  Direct  Numerical  Simulation  approach. 

Probably  the  most  directly  relevant  to  Air  Force  applications  aspect  of  the  present  research  has  been 
the  first-ever  linear  instability  analyses  of  full-bay  cavity  models.  The  key  result  obtained  was  the  (not 
altogether  unexpected)  demonstration  that  the  cavity  geometry  exerts  a  crucial  influence  on  the  principal 
instability  characteristics  of  cavity  eigendisturbances,  namely  their  frequency,  amplification  rate  and  spatial 
structure  of  the  modes’  amplitude  functions. 

In  this  manner,  at  the  conclusion  of  the  Grant,  we  find  ourselves  a  good  deal  closer  to  the  stated 
ultimate  objective  of  our  research,  namely  the  proposition  of  theoretically-founded  flow-control  strategies  for 
realistic  geometries  and  flow  conditions  relevant  to  Air  Force  applications.  The  present  work  has  been  the 
first  to  provide  essential  building  blocks  of  such  a  strategy,  which  permit  relaxing  past  restrictive  assumptions 
on  the  form  of  the  basic  cavity  flow,  and  could  be  integrated  in  flow-control  strategies  aiming  at  addressing 
real  Air  Force  needs. 
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II.  BACKGROUND  AND  MOTIVATION 

Our  concern  in  the  framework  of  this  Grant  has  been  with  hydrodynamic  and  aeroacoustic  instabilities 
in  compressible  open  cavity  flow.  Most  theoretical  work  on  this  problem  to-date  has  focused  on  a  distinct 
feature  of  this  flow,  namely  flow/acoustic  resonance  (Powell  1961;  Rossiter  1964)  through  a  feedback  loop 
involving  (i)  excitation  of  the  unstable  shear  layer  emanating  from  the  upstream  corner  of  the  cavity  through 
a  Kelvin-Helmholtz  (KH)  instability,  (ii)  generation  of  an  unsteady  irrotational  field  by  interaction  of  the 
shear  layer  with  a  solid  boundary,  typically  the  downstream  corner  of  the  cavity  and  (iii)  the  upstream 
influence  of  the  irrotational  field  which  provides  for  the  further  excitation  of  the  instabilities  in  the  shear 
layer,  especially  near  the  upstream  edge.  For  incompressible  flow,  the  upstream  influence  is  instantaneous, 
while  for  compressible  flow  there  is  an  acoustic  delay.  Resonance  occurs  when  the  phase  change  of  a  distur¬ 
bance,  at  a  given  frequency,  leads  to  constructive  reinforcement  and,  ultimately,  saturation.  Because  of  the 
coupling  with  flow  instability,  cavity  resonant  frequencies  and  amplitudes  may  depend  on  the  flow  speed,  U, 
the  boundary  layer  (momentum)  thickness,  0,  just  upstream  of  the  cavity,  ambient  density,  g^,  viscosity,  ju^, 
and  sound  speed,  a*,.  For  a  two-dimensional  rectangular  cavity  with  length,  L,  and  depth,  D,  there  are  thus 
four  dimensionless  parameters  governing  the  flow  L/D,  L/0,  Re0  =  g^UO/ju^ ,  Ma  =  U/  a^  as  well  as  ad¬ 
ditional  parameters  and  profiles  that  would  specify  the  state  of  the  boundary  layer  upstream  of  the  cavity  if  it 
were  turbulent.  The  instability  mechanism  proposed  by  Rossiter  (i.e.  a  convectively  unstable  nearly-parallel, 
two-dimensional  shear  layer  that  is  forced  by  the  large  acoustic  field  produced  near  the  downstream  in¬ 
teraction)  is  but  one  of  several  modes  of  instability  that  have  been  observed  in  flows  over  cavities. 

Factors  that  determine  what  instability  mechanism  will  be  observed  include  all  the  relevant  parameters 
discussed  above.  Even  for  shallow  cavities  (L/D  >  1),  Gharib  and  Roshko  (1987)  observed  in  their 
incompressible  experiments  on  an  axisymmetric  cavity  that,  as  the  length  of  the  cavity  (relative  to  the 
upstream  boundary  layer  thickness)  was  increased,  there  was  a  substantial  change  in  the  behavior  of  the  cavity 
oscillations,  denoted  as  “wake  mode”  by  these  authors.  Under  these  conditions,  the  flow  was  characterized  by 
a  large  scale  (dimensions  of  the  cavity  depth)  shedding  from  the  cavity  leading  edge.  In  two-dimensional 
Direct  Numerical  Simulation  (DNS),  Rowley,  Colonius  and  Basu  (2002)  have  observed  a  very  similar 
transition.  In  both  experiment  and  computation,  the  boundary  layer  upstream  of  the  cavity  was  laminar.  The 
wake  mode  transition  has  also  been  observed  in  two-dimensional  Reynolds  averaged  Navier-Stokes  cal¬ 
culations  at  higher  Reynolds  numbers  (Fuglsang  and  Cain  1992;  Cain,  Rubio,  Bortz,  Banks  and  Smith  2000; 
Shieh  and  Morris  2001).  Rowley  et  al,  (2002)  have  performed  an  extensive  mapping  of  the  transition  to  wake 
mode  for  two-dimensional  flows  in  parameter  space.  A  key  discovery  was  that  the  frequency  of  oscillation  in 
wake  mode  was  found  to  be  nearly  constant  with  Mach  number  (from  0.4  to  0.8)  and  it  would  appear  that  the 
instability  is  in  this  regime  is  purely  hydrodynamic  in  nature. 

Recently  Alvarez,  Kerschen  and  Tumin  (2004)  have  made  progress  in  developing  improved  linear 
models  of  the  scattering/receptivity,  and  acoustic  propagation  aspects  of  the  Rossiter  modes  and  are  able  to 
eliminate  many  of  the  empirical  parameters  in  the  original  Rossiter  model.  However,  when  significant 
interactions  between  the  shear  layer  and  the  flow  within  the  cavity  exist,  the  flow  may  be  found  in  a  state 
completely  different  from  that  of  cavity  resonance  that  can  be  described  with  the  aid  of  the  Rossiter  model 
(or  any  of  its  improved  versions),  namely  the  wake  mode  of  Gharib  and  Roshko  (1987).  The  wake  mode 
appears  to  be  a  global  instability  more  akin  to  vortex  shedding  from  bluff  bodies;  unlike  the  Rossiter  mode 
the  frequency  of  oscillation  is  independent  of  Mach  number,  indicating  that  acoustic  feedback  does  not  play  a 
role  in  the  instability.  However,  wake  mode  is  not  commonly  observed  in  shallow  cavity  flows  at  moderate 
subsonic  to  supersonic  Mach  number.  Further,  it  is  clear  that  three  dimensionality  plays  a  key  role.  Work  by 
Shieh  and  Morris  (2000;  2001)  has  shown  that  two  dimensional  cavities  that  are  oscillating  in  wake  mode 
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return  to  the  Rossiter  mode  when  the  incoming  boundary  layer  has  three  dimensional  disturbances;  it  is  now 
generally  accepted  that  the  development  of  three-dimensionality  in  the  resonant  oscillations  plays  a  role  in 
suppressing  the  wake  mode  instability  (Gloerfelt,  Bogey,  Bailly,  and  Juve  2002).  On  the  other  hand,  it  is 
remarkable,  indeed,  that  the  wake-mode  oscillation  seen  in  the  two-dimensional  simulations,  is  very  similar  to 
that  which  exists  for  deeper  cavities  that  exhibit  higher  levels  of  acoustic  forcing,  such  as  the  oscillations  seen 
in  gas  flows  in  pipes  with  closed  side  branches  (Kriesels,  Peters,  Hirshberg,  Wijnands,  Iafrati,  Riccardi,  Piva, 
and  Bruggeman  1995).  These  are  high-Reynolds  number  three-dimensional  turbulent  flows  developing  in  a 
three-dimensional  geometry.  While  such  modes  have  apparently  not  been  observed  in  shallow  cavities  at  high 
Reynolds  numbers,  it  is  at  present  unclear  under  what  flow  conditions  the  transition  should  occur.  Between 
the  two  extremes  (Rossiter  and  wake  mode)  there  is  a  large  number  of  experimental  observations  that  involve 
significant  interactions  between  the  shear  layer  dynamics  and  the  flow  within  the  cavity  (e.g.  Gloerfelt  et  al. 
2002;  Forestier,  Jacquin  and  Geffroy  2003). 

However,  there  is  very  little  theory  available  with  which  to  understand  all  instability  modes  of  the 
cavity,  their  interaction  and  the  potential  impact  of  such  interaction  on  quantities  of  interest  such  as 
frequency  and  amplitude  of  oscillations,  interactions  of  multiple  tones,  loads  on  internal  cavity  surfaces,  e.t.c. 
The  combined  efforts  in  the  framework  of  the  present  and  the  associated  Caltech  Grants  have  aimed  at  filling 
this  gap  by  calculating  and  understanding  the  instabilities  and  resonance  of  open  cavities  over  the  entire 
parameter  space,  including  those  regimes  where  significant  shear  layer  and  cavity  interactions  occur.  A  key 
element  differentiating  the  current  from  previous  theoretical  efforts  is  the  relaxation  of  the  invariably  used,  in 
past  stability  analyses,  parallel  flow  assumption.  In  the  framework  of  the  present  effort,  the  eigenspectrum  of 
the  entire  cavity  steady  flowfield  has  to  be  computed.  This  flowfield  is  inhomogeneous  in  at  least  two 
(streamwise  and  depth)  spatial  directions,  making  the  open  cavity  an  ideal  candidate  for  the  application  of  the 
BiGlobal  linear  theory  concept,  as  discussed  by  Theofilis  (2003). 

Encouraging  results  in  related  flows  have  recently  appeared  in  the  literature.  Work  by  Barkley,  Gomes 
and  Henderson  (2002)  on  the  global  linear  instability  of  incompressible  flow  in  the  backward-facing  step 
geometry,  which  can  be  related  with  that  of  a  wide  open  cavity  at  appropriate  parameter  ranges,  has 
discovered  the  existence  of  amplified  large-scale  small-amplitude  spanwise-periodic  vortical  structures,  which 
are  distinct  from  the  small-scale  KH-type  of  instability  of  the  shear  layer  emanating  from  the  lip  of  the  step. 
This  global  instability  is  attributed  to  self- excitation  and  amplification  of  an  eigenmode  of  the  closed 
recirculation  bubble  formed  between  the  lip  and  the  downstream  floor  of  the  backward-facing  step.  An 
analogous  mechanism  has  been  proposed  for  the  forward- facing  step  by  S  tuer,  Gyr  and  Kinzelbach  (2000). 
On  the  other  hand,  two  different  models  of  closed  recirculation  bubbles  were  independently  shown  by 
Hammond  and  Redekopp  (1998)  and  Theofilis,  Hein  and  Dallmann  (2000)  to  support  global  linear  instability, 
qualitatively  analogous  to  that  of  flow  in  the  backward- facing  step.  In  particular,  the  latter  work  was  the  first 
to  demonstrate  the  potential  of  a  closed  recirculation  zone  (a  laminar  separation  “bubble”)  to  become  self- 
excited  on  account  of  a  global  instability,  besides  strongly  amplifying  incoming  KH-type  of  disturbances. 
Recent  work  on  this  topic  by  Simens,  Gonzalez,  Theofilis  and  Gomez-Bianco  (2005)  has  established  the 
capacity  of  BiGlobal  linear  theory  to  describe  both  classic  (KH)  and  BiGlobal  eigenmodes  within  a  unified 
framework.  In  cavity  flows,  another  interesting  observation  has  been  made  in  the  experimental  and  numerical 
work  of  Jackson,  Hillier  and  Soltani  (2001)  in  two-dimensional  hypersonic  closed  cavity  flow,  where  Taylor- 
Gortler-like  vortices  forming  in  the  floor  of  the  cavity  appear  to  be  responsible  for  unsteadiness  and  three- 
dimensionality  in  the  flow  at  low  aspect  ratios.  As  the  aspect  ratio  increases  the  structure  spreads  and  covers 
the  entire  cavity,  while  KH  instability  of  the  shear-layer  only  influences  the  laminar-turbulent  transition 
process  at  very  large  cavity  aspect  ratios.  The  global  linear  instability  mechanism  of  flow  in  the  related  lid- 
driven  cavity  model  of  a  closed  cavity  flow  is  now  well-understood  from  a  theoretical  point  of  view  (Theofilis 
2000;  Theofilis,  Duck  and  Owen  2004)  and  also  has  been  verified  experimentally  (Benson  and  Aidun  1992). 
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Evidence  on  the  existence  of  global  modes  in  open  cavities,  prior  to  the  first  demonstration  of  their 
existence  in  the  framework  of  the  research  reported  herein,  could  be  found  in  the  literature.  In  several 
parameter  regimes  of  an  open  cavity,  modes  of  oscillation  that  exhibit  features  of  both  Rossiter  mode  and  the 
wake  mode  have  been  observed.  For  a  cavity  with  Ma  =  0.8  and  L/D  =  0.42  (turbulent  upstream  boundary 
layer)  Schlieren  photographs  by  Forestier  et  al.  (2003)  show  a  flow  structure  that  is  clearly  not  describable  in 
terms  of  the  KH  instability  of  a  nearly  parallel  shear  layer.  Vortex  roll-up  near  the  leading  edge  resembles 
more  closely  that  in  the  early  stages  of  the  wake  mode  cycle,  as  vortical  structures  are  ejected  into  the  free 
stream  to  a  distance  above  the  cavity  which  is  a  significant  fraction  of  the  cavity  length.  DNS  of  two 
dimensional  cavities  shows  that  in  an  intermediate  parameter  regime,  the  flow  may  exhibit  wake  mode 
intermittently,  with  relatively  quiet  periods  dominated  by  shear  layer  instabilities  in  between. 

The  present  theoretical  effort  rests  on  questioning  a  fundamental  assumption  of  the  vast  majority  of 
past  cavity  instability  research.  In  analyzing  the  behavior  of  the  shear  layer,  most  investigators  to-date  have 
implicitly  assumed  that  the  shear  layer  behavior  can  be  described  in  isolation  from  the  geometric  features,  i.e. 
as  if  it  were  a  free  shear  layer.  Nonparallel  effects,  three-dimensionality,  and  the  coupling  of  the  flow  inside 
the  cavity  are  of  course  neglected  in  such  analyses.  An  alternative  analysis  of  the  global  instability  modes, 
which  requires  spatial  DNS  as  input,  has  been  developed  for  the  cavity  by  Theofilis  (2000;  2003).  The  tool  uti¬ 
lized  to  accomplish  this  task  is  the  appropriate  (in  view  of  the  inhomogeneous  geometry)  BiGlobal  linear 
instability  theory.  In  order  to  arrive  at  a  unified  stability  theory  of  cavity  flows,  both  the  compressible 
residuals  algorithm  and  direct  solutions  of  the  partial  derivative  EVP  governing  BiGlobal  linear  instability  of 
nonparallel  two-dimensional  steady  incompressible  and  compressible  flows,  have  been  utilized.  As  with  any 
instability  analysis,  the  first  step  has  been  calculation  of  the  steady  basic  flow  which,  in  the  context  of 
BiGlobal  analysis,  is  obtained  by  two-dimensional  DNS  run  to  steady  state  (if  one  exists  —  more  on  this  issue 
in  Section  III. 2).  In  a  second  step,  two-  and  three-dimensional  instabilities  may  be  calculated  by  numerical 
solution  of  the  appropriate  (partial-differential-equation-based)  eigenvalue  problem  (EVP).  Additionally,  in 
view  of  the  multi-parametric  nature  of  the  problem  indicated  above,  efficient  tools  that  minimize  the 
computing  work  within  the  two-step  approach  discussed  are  needed.  In  the  framework  of  the  present  Grant, 
we  have  extended  the  (EOARD-supported:  F61775-99-W-E090  http:/ /handle.dtic.mil/100.2/ADA379181) 
incompressible  “residuals”  algorithm  to  compressible  flow;  the  compressible  residuals  algorithm  then  formed 
one  of  the  pillars  of  the  Caltech  investigations. 

A  major  thrust  of  the  present  Grant,  not  contemplated  in  the  original  proposal,  but  of  sufficient 
interest  in  order  to  be  pursued  for  over  one-third  of  the  duration  of  the  Grant,  has  been  the  (worldwide  first) 
investigation  of  the  effect  on  instability  characteristics  (frequency  and  spatial  distribution  of  the  amplitude 
functions  of  the  eigenmodes)  of  an  (empty)  open  cavity,  that  objects  placed  inside  the  cavity  may  have.  The 
technology  developed  herein  has  intended  to,  and  to  a  large  degree  succeeded  in  merging  the  gap  between  the 
academic  research  capabilities  existing  prior  to  the  present  effort  and  the  actual  warfighter  needs,  in  terms  of 
direct  prediction  of  eigenfrequencies  (tones)  of  complex  full-bay  cavity  configurations.  At  the  conclusion  of 
the  Grant,  we  are  in  a  position  to  perform  BiGlobal  instability  analyses  of  relatively  complex  full-bay  model 
configurations,  under  the  sole  assumption  of  a  steady  laminar  (non-parallel)  two-dimensional  basic  state. 
Analyses  of  time-periodic  laminar  or  Reynolds-averaged  turbulent  basic  flows  are  straightforward  extensions 
of  the  present  work,  which  would  make  use  of  the  analysis  tools  developed  herein.  In  what  follows,  the 
theoretical  approaches  used  are  briefly  exposed  in  Section  III,  making  reference  to  material  published  in  the 
course  of  the  present  investigations  for  further  details.  In  Section  IV  results  obtained  during  the  execution  of 
the  Grant  are  outlined,  while  emphasis  is  placed  on  the  most  novel  development  aspects  of  the  present 
research,  namely  the  spectral  multi-domain  algorithm  which  permits  study  of  instability  in  complex  cavity 
configurations.  Our  conclusions  and  projected  outlook  on  future  research  is  offered  in  Section  V. 
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III.  THEORETICAL  APPROACHES 

1.  The  compressible  BiGlobal  eigenvalue  problem  (EVP) 

In  the  course  of  the  Grant  concluded,  BiGlobal  instability  theory  was  applied  for  the  first  time  to  open 
cavity  flows  in  both  the  incompressible  and  the  compressible  flow  regimes.  A  review  of  this  analysis 
methodology  and  its  application  to  a  wide  spectrum  of  flows  has  been  presented  by  Theofilis  (2003).  The 
analysis  proceeds  with  the  compressible  equations  of  motion,  from  which  the  incompressible  BiGlobal  EVP 
may  also  be  recovered.  The  equations,  written  for  an  arbitrary  coordinate  system,  are 

dp 


dt 


+  V  •  (pu)  =  0, 


d(pu ) 
dt 
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where  shorthand  notations,  respectively,  for  the  viscous  stress  tensor  a, 


o  =  p  |  (Vu  +  VuT)  -  ^(V  •  u)  I 


and  the  viscous  dissipation  function  ®, 


<I>  =  I :  V  u  —  V  ur  I  :  a 


have  been  used. 


Linear  stability  analysis  in  the  BiGlobal  framework  involves  the  substitution  of  a  decomposition  of  any 
of  the  independent  flow  variables,  e.g.  the  three  velocity  components,  temperature  and  pressure  q = 
(u,  v,  w,  0,  p)T  into  the  equations  of  motion.  All  quantities  are  considered  to  be  composed  of  an  0(1)  steady 
two-dimensional  basic  state  and  small-amplitude  0(e)  unsteady  three-dimensional  perturbations,  according  to 
the  BiGlobal  Ansatz  shown  in  Table  1.  In  this  table,  single  primes  refer  to  time-  or  space-periodic  quantities, 
and  double  primes  refer  to  spatial  coordinates  with  mild  variation  in  comparison  with  that  on  the  non-primed 
spatial  coordinates.  Note  that  the  BiGlobal  decomposition  shown  in  Table  1  is  applicable  to  both 
incompressible  and  compressible  flows.  In  the  compressible  flow  at  hand,  upon  substitution  into  the  general 
equations  of  motion  of  the  BiGlobal  instability  Ansatz,  use  of  a  perfect-gas  constitutive  law  and  linearization, 
the  compressible  non-symmetric  generalized  linear  BiGlobal  EVP 

£  q  =  Q"R-  q, 

for  the  determination  of  the  complex  eigenvalue  Q  =  Qr  +  i  Qi  results.  In  the  temporal  framework  used,  Qr 
represents  a  frequency  and  Qi  the  amplification/ damping  rate  of  the  disturbance  sought. 
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Table  1:  Linear  Stability  Theory  Decompositions. 
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Explicitly,  the  EVP  reads 
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where  the  unknown  vector  in  this  formulation  comprises  disturbance  velocity  components  and  disturbance 
temperature,  respectively  corresponding  to  solution  of  the  momentum  and  energy  equations,  and  disturbance 
pressure,  corresponding  to  the  equation  of  continuity.  The  lack  of  boundary  conditions  for  the  latter  quantity 
requires  the  use  of  staggered  grids  and  a  natural  choice  ensuring  high  accuracy  on  the  minimum  number  of 
nodes  possible  is  spectral  collocation  (Macaraeg,  Streett  and  Hussaini  1988).  Spectral  collocation  based  on 
two-dimensional  mapped  Chebyshev  Gauss-Lobatto  (CGL)  and  Chebyshev  Gauss  (CG)  grids  is  used  for 
spatial  discretization  of  the  EVP,  as  shown  in  Figure  1. 


13 


FINAL  REPORT 

Grant  FA8655-03-1-3059  (Theofilis)  -  UA  Unified  View  of  Global  Instability  of  Compressible  Flow  over  Open  Cavities” 


Figure  1:  Upper:  The  Chebvshev  Gauss-Lobatto  (CGL)  and  the  Chebvshev  Gauss  (CG)  points,  respectively  denoted  by  full  and 
open  circles.  Lower:  the  resulting  discretization  on  a  two-dimensional  CGL  and  a  staggered  two-dimensional  CG  grid. 
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The  CGL  grid  is  used  in  order  to  collocate  the  disturbance  momentum  and  energy  equations,  while  the  CG 
grid  is  used  in  order  to  collocate  the  disturbance  continuity.  The  operators 


T  —  TG  'T  —  tgl 
x  —  IGL^  J  ~  1G 


denote  interpolation  between  the  Chebyshev  Gauss-Lobatto  and  Chebyshev  Gauss  grids  and  are  used  in 
order  to  transfer  information  between  the  two  grids. 


The  linear  operators  are  defined  in  the  Appendix.  Solution  methods  for  the  eigenvalue  problem  are 
discussed  in  detail  by  Theofilis  (2003);  a  pictorial  representation  of  the  different  paths  followed  is  shown  in 
Figure  2.  Unless  validation  studies  are  performed,  for  which  the  direct  QZ  algorithm  is  used  in  order  to 
recover  the  full  eigenspectrum  at  a  cost  scaling  with  the  cube  of  the  total  number  of  nodes  used  in  the 
discretization,  the  bulk  of  the  analyses  is  performed  using  Krylov  subspace  iteration  and  the  (implicitly 
restarted)  Arnoldi  method  (IRAM).  However,  it  should  be  noted  that  the  matrix 

e  =  % 


discretizing  the  EVP  may  be  formed  and  kept  in  memory,  in  single-domain  computations.  This  is  not 
anymore  the  case  in  the  multidomain  solutions  of  the  EVP,  as  will  be  discussed  in  Section  IV.6. 
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Figure  2:  Schematic  representation  of  the  possible  numerical  procedures  for  the  solution  of  the  BiGlobal  EVP. 


2.  The  compressible  residuals  algorithm 

The  algorithm  rests  on  the  identification  of  small-amplitude  residuals  in  the  (unsteady)  DNS  calculation 
for  the  recovery  of  a  steady  state  as  the  least  damped  global  two-dimensional  eigenmodes  of  this  underlying 
steady  flow  (http :  /  /  handle .  dtic .  mil  7100.2/  AD  A3 79181).  Concretely,  referring  to  Table  1,  the  pertinent 
decomposition  is  that  of  BiGlobal  theory  in  the  limit  of  P=0.  It  is  then  possible  to  utilize  this  information  in 
order  to  recover  the  converged  steady-state  solution  as  well  as  the  amplitude  functions  of  the  least-damped 
eigenmode  from  transient  data,  without  having  to  pursue  the  time  integration  of  the  equations  of  motion 
until  convergence  of  the  flow  in  time  is  obtained. 


Consider  a  transient  solution  q  of  the  equations  of  motion  near  convergence  to  the  respective  steady 
state,  denoted  by  the  barred  quantity,  as  containing  a  small-amplitude  two-dimensional  disturbance  according 
to: 

q(x,  y,  t)  =  qfz,  y)  +  £  q2D(:r,  yJ), 

where 

q2n(x,  y,t)  -  q2D(T  ;y)  Q<Tt 


Combining  this  expression  with  the  BiGlobal  theory  Ansatz  shown  in  Table  1,  one  obtains 


q(.r,  y ,  £)  —  q(x,  y)  +  £  q,  cos —  qi  sino>r£ 


where  G  may  be  defined  from 

a  =  ln[qf/qf-^]/A£  «  dln^j/dt. 
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with 

<lf  =  |q(*o,lM)  -  q(*of  yo)|. 

The  calculation  of  the  basic  state  and  the  amplitude  functions  from  transient  data  follows  in  two  stages.  First, 
elementary  signal  analysis  techniques  deliver  the  results  for  the  circular  frequency  <Dr.  The  circular  frequency  is 
calculated  from  the  period  of  oscillations  in  the  transient  time-signal  of  q  which,  in  turn,  is  identified  by  the 
maxima  in  the  signal.  Second,  once  both  <Dr  and  G  have  converged  in  time,  the  BiGlobal  Ansatz  may  be  used 
to  calculate  both  the  steady  basic  state  and  the  amplitude  functions  of  the  least-damped  (P=0)-eigenmode. 


Two  independent  paths  may  be  used  in  order  to  calculate  G,  depending  on  whether  the  transient  signal 
is  oscillatory  or  not.  In  the  first  case,  the  equation 


W  +  la  +u')m~-aW 


=  (J. 


may  be  used,  with  0Dr.  and  the  time-derivatives  known.  This  expression  may  be  evaluated  at  those  times  that 
the  first  time-derivative  of  the  signal  is  zero  in  the  course  of  the  time-integration,  i.e.  at  the  same  times  that  cor 
is  calculated.  At  these  times  the  magnitude  of  G  is  given  by 


1  (33q  j&P) 

2  (£2q jdf2) 


(dqfdt)= 0 


On  the  other  hand,  in  case  G0r  =  0,  a  monotonic  dependence  of  the  first  time  derivative  of  the  signal  on 
time  is  typically  observed  from  the  beginning  of  the  calculation  until  convergence,  with  the  time-derivative  of 
the  signal  being  zero  possibly  during  early  transients  and  certainly  at  convergence.  Avoiding  calculations 
during  the  early  transient  stages,  the  magnitude  of  G  may  be  calculated  using 


(gy  ot2) 

(dq/dt) 


With  <Dr  and  G  known,  the  BiGlobal  Ansatz  may  be  written  for  the  case  (P=0)  as  a  (complex)  linear 
system  of  three  equations  at  three  times,  q  =  t,  t2=t  +  At,  and  t3  =  t  +  2At,  for  three  unknowns,  namely  the 
steady  basic  state,  and  the  real  and  imaginary  parts  of  the  amplitude  function  of  the  least-damped  eigenmode, 
with  the  transient  solution  known  at  these  three  times.  Simple  algebra  then  delivers  the  desired  converged 
steady-state  solution  as 

qi  e2crif  —  2  c[2  eaAt  cos  At  4-  qa 
I  e2crAf  _  2  effAf  cos  u)r At  +  1 


With  the  steady-state  solution  known,  the  amplitude  function  of  the  least-damped  eigenmode  may  be 
recovered  from 
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A 

ft  - 

si(<i2  -  q)  -  ss(qi  -  q) 

%  — 

C2$l  -  ci  $2 

A 

qi  = 

ci(q-2  -  q)  -  c2(qi  -  q) 

C2'sl  —  Cl  $2 

where 

Cl  = 

expat].  cos 

II 

exp  at 2  cos  ojTt2 

Si  — 

exp  ati  sin  topi 

^2  = 

exp  oU  sin  ■ 

The  accuracy  of  the  results  delivered  by  the  residuals  algorithm  clearly  depends  on  that  by  which  the 
circular  frequency  and  the  damping  rate  are  provided,  which,  in  turn  depends  on  the  accuracy  by  which  the 
first  three  time-derivatives  of  the  transient  solution  are  calculated.  The  latter  is  a  function  of  the  time-step  in 
the  calculation  and  the  number  of  fields  stored,  in  order  for  backward  differentiation  formulae  to  be  applied. 
Since  the  time-step  is  controlled  by  CFL  considerations,  it  is  advisable  to  store  a  reasonably  high  number  of 
fields  in  order  for  high  accuracy  of  cor  and  G  to  be  obtained. 


3.  A  spectral  multi-domain  algorithm  for  the  BiGlobal  EVP  in  complex  cavity  configurations 

As  mentioned,  the  residuals  algorithm  was  made  available  to  the  collaborative  Caltech  research  and  has 
been  used  a  good  part  of  their  investigations.  The  present  effort  was  then  diverted  toward  direct  solution  of 
the  BiGlobal  EVP,  for  both  the  (empty)  open  cavity  and  full-bay  models.  In  the  latter  case,  a  new  multi- 
domain  algorithm  was  designed,  validated  and  applied  to  various  open  cavity  configurations,  as  discussed  in 
detail  in  Section  IV.  Here  the  essential  aspects  of  spatial  discretization  are  exposed.  For  ease  of  presentation, 
the  incompressible  flow  limit  is  considered,  where  the  BiGlobal  EVP  reads 

(LHS)  U  =  A  (RHS)  U , 

with 


U 

L2-{D,(u )) 

-dm, 

0 

-V 

(I 

0 

0 

(T 
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Here 

L1  =  (d)  +  D1y-/31)-uDy-vDy 

barred  and  hatted  quantities  denote  basic  and  disturbance  flow  quantities,  respectively,  and  D ■  denotes 
differentiation  w.r.t.  the  i-th  spatial  dimension.  Multi-domain  resolution,  schematically  depicted  in  Figure  3, 
involves  an  extended  linear  system  as  follows: 


Figure  3:  Multi-domain  spatial  discretization  of  a  square  domain  by  two  rectangular  sub-domains,  denoted  “green”  and  “red”. 


where  LHSred  and  LHSgreen  are  linear  operators  of  the  same  structure  as  that  describing  the  single¬ 
domain  BiGlobal  EVP  solution  and  Cl  denotes  compatibility  conditions  at  interfaces,  necessary  in  order  to 
couple  the  two-dimensional  amplitude  functions  of  the  eigenmodes  in  the  different  domains.  From  the 
layout  of  the  multi-domain  matrix  LHS  it  becomes  evident  that  the  direct  methods  used  for  the  residual 
matrix  operations  (matrix  factorization,  matrix-vector  multiplication)  in  the  work  of  Theofilis  et  al.  (2004) 
become  increasingly  inappropriate  for  the  solution  of  the  spectral  multi-domain  problem  as  the  number  of 
sub-domains  increases.  The  reason  is  that  the  leading  dimension  of  the  multi-domain  matrix  LHS  is  a 
multiple  of  that  of  the  single-domain  matrices  (already  very  large)  LHSred and  LHS green.  Consequently,  particular 
efforts  were  devoted  toward  efficient  solutions  of  the  multi-domain  EVP,  the  details  of  which  are  provided  in 
Section  IV. 
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IV.  TECHNICAL  RESULTS 

1.  Computing  cluster  installation  and  performance 

In  preparation  for  the  execution  of  the  project,  computing  hardware  has  been  acquired,  according  to 
the  provisions  of  the  Grant.  Such  hardware  corresponds  to  an  8-node  cluster  (16  processors/ 16Gbytes),  as 
originally  planned.  However,  the  drop  in  prices  between  the  time  at  which  the  Grant  has  been  approved  and 
that  of  the  hardware  acquisition  has  resulted  in  the  ability  to  purchase  a  fast  (Myrinet)  inter-nodal  connection 
hardware,  to  ensure  rapid  communication  of  data  in  the  invariably  parallel  computations  within  the  project,  in 
both  the  direct  numerical  simulations  -  DNS  and  the  BiGlobal  stability  analysis  aspects  of  the  work. 
Furthermore,  rather  than  investing  in  a  stand-alone  cluster  facility,  the  hardware  configuration  acquired  was 
chosen  to  be  compatible  with  a  newly-built  computing  cluster  facility  at  the  Universidad  Politecnica  de 
Madrid  (UPM),  the  academic  affiliation  of  Theofilis.  This  cluster  currently  comprises  a  total  of  64  nodes  (128 
processors/256  Gbytes)  and  there  exist  15  Tbytes  of  storage  space.  The  cluster  runs  under  Linux, 
parallelization  is  done  using  MPICH  and  PBS  is  used  as  job-scheduling  software.  The  decision  to  invest  in 
compatible  hardware  was  taken  with  the  objective  in  mind  to  maximize  the  computing  power  available  to  nu- 
modelling,  beyond  the  8-node  configuration  pertaining  to  the  company.  The  modus  operandi  during  the 
execution  of  the  Grant  has  been  that,  when  the  8-nodes  pertaining  to  nu-modelling  were  not  needed  in  their 
entirety  (e.g.  during  code-development  phases),  they  have  been  used  by  other  UPM  members.  In  exchange, 
when  production-runs  associated  with  the  Grant  have  been  necessary,  nu-modelling  had  access  to  a  computing 
facility  that  is  substantially  larger  than  that  owned  by  the  company.  The  decision  to  spend  a  substantial 
portion  of  the  funds  set  aside  for  computing  hardware  in  inter-nodal  connection  hardware,  as  opposed  to  a 
larger  number  of  processors,  has  been  taken  on  grounds  of  hardware  performance,  as  explained  next. 

Commensurate  with  the  first  objective  of  the  first  phase  of  the  Grant,  the  (MPICH-)  parallelized  2d- 
DNS  code  of  Caltech  has  been  ported  on  this,  as  well  as  predecessors  of  the  aforementioned  cluster.  Speed¬ 
up  factors  have  been  realized  when  using  the  (’’new”)  cluster,  compared  with  performance  on  a  machine 
available  at  UPM  prior  to  this  cluster,  designated  as  “old”  in  what  follows.  Speed-up  has  been  realized 
commensurate  with  two  factors,  the  processor  and  the  communication  speed.  Table  2  shows  a  comparative 
study  of  performance  of  the  same  version  of  the  DNS  code.  Shown  are  the  processor  type,  Proc,  the  time 
spent  while  computing,  tCPU,  and  the  time  spent  in  data-exchange  (communication),  tCOM,  both  in  seconds  per 
step  in  the  code,  the  actual  (as  opposed  to  theoretical)  sum  of  these  two  times,  tSTEP,  in  seconds  and  the  ratio 
tCOM  /  tsTEP?  indicating  the  percentage  of  time  spent  communicating  data,  as  opposed  to  computing.  In  order 
to  facilitate  comparisons,  the  same  number  of  20  processors  has  been  used  on  both  the  “old”  and  the  “new” 
clusters.  Both  machines  use  fast  inter-node  connection  (Scali  in  the  “old”  and  Myrinet  in  the  “new”  cluster). 
For  this  moderate  computation  the  speed-up  realized  on  the  “new”  machine  is  due  to  both  its  faster 
processor  type  and  its  faster  inter-nodal  connection  hardware. 

Table  2:  Comparative  study  of  code  performance  on  the  two  clusters. 


Cluster 

Proc 

(GHz) 

*CPU 

(s/STEP) 

t-COM 

(s/STEP) 

*STEP 

(s) 

LoM  /  tsTEP 

(%) 

“old” 

P-III  (1.8) 

64.23 

18.60 

88.4 

21 

“new” 

P-IV  (3.06) 

45.97 

8.97 

57.7 

15 
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The  situation  changed,  as  the  computation  became  more  demanding  and  a  larger  number  of  processors 
and  memory  were  necessary.  The  speed-up  factor  in  this  case  turned  out  to  be  determined  by 
communication-,  rather  than  calculation-speed.  Figure  4  shows  the  speed-up  factor,  this  time  measured  in 
terms  of  number  of  steps  realized  by  the  DNS  code  per  second  (S/sec)  as  the  number  of  processors  (np) 
increases.  For  this  comparison  the  “new”  cluster  has  been  used,  communicating  data  either  through  its 
Gigabit  network  (used  for  NFS  operations)  or  through  its  Myrinet  network  (used  for  computations).  It  may 
clearly  be  seen  that  both  solutions  performed  comparably  well  at  low  number  of  processors;  however,  as  np 
increased,  the  Myrinet-based  computation  scaled  perfectly  linearly,  i.e.  the  number  of  time-steps  calculated  increased 
linearly  with  the  number  of  processors  used.  On  the  other  hand,  past  a  threshold,  the  Ethernet-based  solution  reached  a 
plateau  of  number  of  time-steps  computed \  irrespective  of  the  number  of  processors  used.  In  other  words,  an  Ethernet-based 
solution  would  have  become  a  bottleneck  for  large-scale  computations,  of  the  type  actually  performed  later, 
for  both  three-dimensional  multi-domain  BiGlobal  instability  analysis  and  DNS.  With  hindsight,  this  has 
justified  the  decision  to  balance  the  expenditure  of  the  funds  made  available  within  the  scope  of  the  Grant 
between  number  of  nodes  and  (expensive  but  fast)  inter-nodal  connection  Myrinet  hardware. 

Figure  4:  Comparative  study  of  code  performance  on  an  Ethernet-based  (clotted  line)  and  a  Myrinet-based  (solid  line)  cluster. 


A  parallel  development,  unforeseen  at  the  time  of  application  of  the  grant,  has  been  relevant  to  the 
Grant  execution,  namely  the  installation  of  two  large  supercomputer  facilities  in  Spain.  Mare  Nostrum  has  been 
installed  in  the  Barcelona  Supercomputing  Center,  while  its  smaller  counterpart,  CeSViMa  (Centro  de 
Supercomputacion  y  Visualizacion  de  Madrid)  has  been  installed  at  (and  owned  by)  the  Universidad 
Politecnica  de  Madrid;  they  are  both  IBM  Linux  clusters,  respectively  featuring  ~4k  and  ~lk  processors 
interconnected  by  Myrinet  hardware.  The  existence  of  the  first  of  these  machines  has  overlapped  the  longest 
with  the  duration  of  the  Grant  but  has  been  under  testing  throughout  the  entire  same  period,  so  that  no 
access  to  it  was  requested.  The  local  supercomputer  has  been  utilized  toward  the  end  of  the  Grant,  in  order 
to  clarify  parallel  performance  issues  with  respect  to  the  sparse-matrix  solution  algorithms  used  in  the  spectral 
multi-domain  solution  of  the  BiGlobal  EVP.  It  was  shown  that  no  gains  are  to  be  expected  by  performing  the 
analysis  related  with  the  Grant  on  a  supercomputer,  such  that  for  the  bulk  of  the  analyses  reported  below,  the 
numerical  work  has  been  executed  on  the  local  UPM  cluster  described  above. 
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2.  First  application  of  the  compressible  residuals  algorithm  in  compressible  flow 

At  the  beginning  of  the  Grant  execution,  the  connection  with  previously-initiated  global  instability 
analysis  research  was  made  by  applying  for  the  first  time  to  a  compressible  flow  instability  problem  the 
“residuals”  algorithm  devised  by  Theofilis  during  EOARD-sponsored  work  (F61775-99-W-E090 
http:/ /handle. dtic.mil/ 1 00.2 /ADA379 181)  for  incompressible  flows.  This  algorithm  may  be  used  in 
conjunction  with  two-  or  three-dimensional  DNS  in  order  to  recover  the  underlying  instability  characteristics 
of  the  flow  from  transient  DNS  data.  Prior  to  the  Grant,  the  algorithm  had  only  been  applied  to 
incompressible  flows  while  during  the  Grant,  and  in  collaboration  with  the  Caltech  team,  the  algorithm  has 
been  applied  for  the  first  time  to  compressible  DNS  data,  as  discussed  below. 

Technical  details  of  the  residuals  algorithm  were  presented  in  the  Section  III. 2.  In  brief,  the  residuals 
algorithm  uses  simple  algebraic  operations  on  unsteady  DNS  data  and  delivers  both  the  steady-state  and  the 
instability  characteristics  of  any  (one-,  two-,  or  three-dimensional)  flow  from  unsteady  DNS  data,  the  latter 
taken  during  the  early  stages  of  the  integration  time  necessary  to  reach  convergence.  The  alternatively 
available  methodology  in  order  for  instability  results  to  be  obtained  is  to  march  the  unsteady  DNS  until 
convergence  to  a  steady-state  and  subsequently  solve  the  appropriate  corresponding  (one-,  two-  or  three- 
dimensional)  eigenvalue  problem  (EVP).  It  is  known  that  the  cost  of  the  DNS  scales  with  the  powers  6/4 
and  9/4  of  Reynolds  number  in  two-  and  three-spatial  directions,  respectively.  Moreover,  as  discussed  by 
Theofilis  (2003),  the  leading  dimensions  of  the  matrices  involved  in  the  solution  of  the  BiGlobal  and 
TriGlobal  EVP  scale  with  the  product  of  the  degrees  of  freedom  used  to  resolve  two-  or  three  spatial 
directions,  respectively,  such  that  the  residuals  algorithm  permits  recovery  of  the  sought  instability  results  at 
an  insignificantly  small  fraction  of  the  otherwise  necessary  computing  effort.  An  example  of  successful 
application  of  the  algorithm  on  the  classic  lid-driven  cavity  problem  is  shown  in  Table  3  and  the  associated 
Figure  5.  Single  overbar  denotes  time  at  which  the  DNS  solution  converges,  while  double  overbar  denotes 
time  at  which  the  residuals  algorithm  was  applied;  in  all  cases  the  latter  is  of  0(20%)  of  the  former  time, 
indicating  the  substantial  savings  achieved.  The  shorthand  notation 

^3=  l(l-q)/q| 


denotes  the  relative  error  between  quantity  q  evaluated  at  the  aforementioned  early  time  over  the 
respective  converged-in-time  quantity,  the  latter  obtained  at  the  end  of  the  time  integration. 

Table  3:  Recovery  of  the  stream  function  ^  in  the  square  liddriven  cavity  problem,  from  transient  DNS  data  taken  at  time  t=20. 

out  of  a  time-to-convergence  being  t=100.  The  effect  of  resolution  and  timestep  used  in  the  DNS  is  examined  at  Re=100; 

x(v)  =  x  x  !Qy  (Theofilis  and  Colonius  2003). 
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1 
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32  x  32 
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i 
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0.540870 

max(A^) 

5.3{— 8) 

8.8(-7) 

4.0(— 0) 

min(Av) 

3.6{-9) 

7. 9  (—8) 

5.1  (—7) 

max(  Aufi) 

3.4(  —4) 

3.5(— 4) 

1.0(— 3) 

min(A^) 

1G(-G) 

1  2  (—3) 

1_6(— 2) 
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The  most  significant  results  here  are  (i)  the  low  level  of  approximation  error  max(At|>)  and  (ii)  the  ratio 
of  the  (“double-overbar”)  time  at  which  the  transient  solution  is  processed  to  that  (“overbar”)  otherwise 
necessary  for  the  DNS  to  converge.  The  case  Re  =  100  is  typical  of  one  in  which  the  least  stable  eigenmode 
determines  the  transient  behavior  of  the  DNS  throughout  most  of  the  time  integration  process.  With  the 
results  for  the  damping  rate  a  converging  quite  quickly,  the  desired  converged  steady  state  may  be  obtained  at 
a  time  between  a  quarter  at  the  coarsest  and  a  fifth  at  the  finest  resolution  of  the  time  required  by  the  time  - 
marching  algorithm  for  the  residuals  to  be  eliminated.  The  result  for  a  is  only  marginally  affected  by 
resolution  and  time  step;  the  precise  times  at  which  a  converges  are  affected  by  a  small  amount  when  refining 
the  grid,  with  the  finest  resolution  results  converging  earlier.  In  all  cases  application  of  the  residuals  algorithm 
results  in  substantial  savings  compared  with  the  otherwise  necessary  computing  effort. 


Table  4  Parameters  of  compressible  open  cavity  L/D— 2  DNS  runs  performed. 


DNS  Run 

I 

II 

Ree 

56.9 

91.0 

ReD 

1500 

2667 

Ma 

0.2 

0.8 

The  first  application  of  the  residuals  algorithm  to  compressible  flows  was  presented  by  Theofilis  and 
Colonius  (2003).  The  well- validated  high-order  DNS  methodology  of  Caltech  has  been  used  for  the  recovery 
of  the  basic  states  that  were  analyzed,  while  the  residuals  algorithm  was  derived  and  coded  in  a  joint  effort  of 
the  nu  modelling  and  Caltech  teams.  Basic  flow  results  of  two  simulations,  the  parameters  of  which  are 
summarized  in  Table  4,  have  been  isolated  for  some  discussion  here,  owing  to  their  representative  nature. 
The  time-dependence  of  the  normal  velocity  component  at  a  single  location  inside  the  cavity,  shown  in  Figure 
6,  is  characteristic  of  the  signal  returned  by  the  DNS  I  for  any  flow  quantity.  Such  a  signal  is  composed  of 
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two  damping  BiGlobal  eigenmodes,  a  traveling  one  visible  in  the  signal  during  t  e  [0,  300]  and  a  stationary 
one,  which  emerges  (although  damped  throughout  the  time-integration)  beyond  t  >  300.  Before  discussing 
the  characteristics  of  the  BiGlobal  instability  eigenmodes,  the  residuals  algorithm  has  first  been  employed  in 
order  to  recover  the  steady  state  towards  which  the  solution  is  converging.  Owing  to  the  presence  of  multiple 
modes  in  the  simulation,  a  relatively  long  time-integration  has  been  necessary  before  the  residuals  algorithm 
could  be  applied.  Nevertheless,  the  comparison  of  the  estimated  and  the  exact  steady  states  revealed 
approximation  errors  of  10"10,  in  line  with  the  analogous  result  in  the  lid-driven  cavity  flow. 


Figure  6:  Time  development  of  the  normal  velocity  component  in  DNS  I 


Ma  =  0.2 


time 

Turning  to  the  spatial  structure  of  the  BiGlobal  eigenmodes  being  damped  in  the  DNS  process,  Figures 
7  and  8  show  the  real  and  imaginary  parts  of  the  complex  amplitude  functions  corresponding  to  the  wall- 
normal  disturbance  velocity  and  the  disturbance  pressure,  as  obtained  at  t  =  100  and  600,  respectively. 
At  t  =  100  the  traveling  eigenmode  is  dominated  by  a  Tollmien-Schlichting  wave  developing  along  the 
downstream  cavity  wall  and  being  damped  in  time.  At  t  =  600  this  mode  has  subsided  and  the  only  remaining 
disturbance  is  a  stationary  mode,  which  needs  to  be  damped  before  the  steady-state  solution  is  obtained. 

Three  properties  of  this  mode,  its  being  stationary,  its  scaling  with  a  geometric  length  scale  of  the 
cavity,  and  its  location  inside  the  open  cavity  point  to  this  BiGlobal  eigendisturbance  being  related  with  the 
wake-mode  of  Gharib  and  Roshko  (1987).  On  the  other  hand,  the  spatial  structure  of  the  amplifying 
eigenmode  manifesting  itself  in  a  linear  manner  at  early  times  during  DNS  II  has  a  distinctly  different  spatial 
structure  of  the  amplitude  functions.  Figure  9  also  shows  the  same  quantities  as  presented  in  Figures  7  and  8; 
in  both  results  one  may  notice  the  following  features.  First,  the  dominant  structure  is  located  at  the  open  end 
of  the  cavity  and  appears  to  be  related  with  instability  in  the  shear-layer,  of  the  class  predicted  by  the  Rossiter 
analysis.  Clearly,  this  disturbance  is  only  part  of  the  amplitude  function,  a  significant  part  of  which  may  be 
identified  as  a  Tollmien-Schlichting  wave  on  the  downstream  wall  of  the  open  cavity.  Furthermore,  acoustic 
pressure  disturbances  may  be  identified  in  the  respective  amplitude  function,  apparently  originating  at  the 
downstream  corner  of  the  cavity  and  propagating  upstream  in  the  field.  The  acoustic  disturbances  are  weak  in 
comparison  with  the  maxima  of  their  respective  hydrodynamic  counterparts.  The  most  interesting  aspect  of 
these  results  is  that  all  instability  features,  which  have  received  attention  independently  in  the  past,  are 
part  of  one  and  the  same  BiGlobal  flow  eigenmode. 
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Figure  7:  Real-  and  imaginary  parts  of  the  wall-normal  complex  disturbance  velocity  component  (upper  row)  and  pressure 

perturbation  (lower  row)  in  DNS  I  at  t=100 


Figure  8:  Same  as  in  figure  7,  at  t— 600 
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Figure  9:  Same  as  in  figure  8.  for  DNS  II 


3.  Classes  of  instability  modes  discovered:  a  unifying  perspective 

The  first  application  of  the  residuals  algorithm  for  the  recovery  of  BiGlobal  instability  results  from 
DNS  data  to  the  compressible  open-cavity  flow  problem  has  yielded  the  spatial  structure  of  the  most 
energetic  members  of  the  BiGlobal  eigenspectrum,  that  are  otherwise  accessible  to  numerically  challenging 
solutions  of  the  partial-derivative  eigenvalue  problem,  or  to  DNS.  The  present,  and  analogous  results  not 
presented  here,  have  permitted  forming  a  unified  point  of  view  of  open  cavity  flow  instabilities.  Known 
shear-layer  (Kelvin-Helmholtz  —  KH),  Tollmien-Schlichting  (TS)  or  newly-discovered  (wake-mode-related) 
hydrodynamic  and  aeroacoustic  instabilities  in  open  cavities  have  been  shown  to  emerge  as  components 
of  the  amplitude  function  of  single  BiGlobal  eigenmodes.  This  re-interpretation  is  broader  than 
previously  used  modeling  approaches  in  an  essential  manner,  namely  in  that  it  permits  addressing  the 
nonparallel  open  cavity  flow  instability  problem  in  its  entirety,  rather  than  focusing  on  one,  or  some  of  its 
features  (such  as  the  shear-layer,  the  downstream  TS  mode  or  the  acoustic  wave  emanating  from  a  corner)  in 
isolation  from  the  rest  of  the  flow.  The  presently  offered  framework  opens  new  modeling  avenues  for  such 
instabilities.  In  line  with  two  of  the  objectives  of  the  Grant,  this  result  has  provided  a  novel  framework  for 
classification  and  a  unified  perspective  of  all  cavity  instabilities.  A  systematic  study  of  this  point  was 
continued  in  the  framework  of  the  related  Caltech  effort.  At  nu  modelings  results  obtained  to  this  point 
underlined  the  necessity  to  develop  a  new  compressible  BiGlobal  eigenvalue  problem  solver,  in  order  to  aid 
the  study  of  open  cavity  flow  disturbances  with  an  essentially  inhomogeneous  nature,  such  as  the  wake-mode 
instability  outlined  above. 
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4.  On  the  role  of  three-dimensionality  in  the  cavity 

Prior  to  embarking  upon  new  developments  towards  satisfying  the  next  objective  of  the  Grant,  namely 
advance  the  current  understanding  of  the  role  of  three-dimensionality  in  cavity  instability,  an  incompressible 
model  of  three-dimensional  cavity  flow  has  been  devised  and  analysed  using  BiGlobal  instability  analysis.  The 
reason  for  introduction  of  this  model  has  been  that,  while  on  occasion  BiGlobal  instability  analysis  delivers 
results  in  excellent  agreement  with  experiment  (e.g.  Theofilis  2000  —  EOARD-sponsored  research  F61775- 
99-W-E090)  there  exist  anomalies  in  the  predictions  of  linear  (both  classical  ID  and  BiGlobal)  theory.  Most 
notable  of  these  anomalies,  explainable  within  the  framework  of  transient  growth  theory  (Schmid  and 
Henningson  (2001),  is  the  failure  of  ID  linear  theory  to  predict  instability  in  Hagen-Poiseuille  (pipe)  flow, 
aggravated  by  the  erroneous  predictions  of  BiGlobal  theory  of  stability  of  pressure-gradient  driven  flow  in  a 
square  duct  (Tatsumi  and  Yoshimura  1990)  and  wall-bounded  Couette  flow.  It  thus  becomes  of  interest, 
before  investing  efforts  in  the  development  of  a  new  compressible  BiGlobal  EVP  solver  in  an  open  cavity 
configuration,  to  provide  (even  partial)  answers  to  the  question  whether  application  of  Bi-/TriGlobal  linear 
theory  to  the  open  cavity  is  expected  to  deliver  realistic  results.  Resolution  of  this  issue  has  been  particularly 
pressing  at  that  stage  of  the  Grant,  since  the  conclusions  reached  on  the  aforementioned  question  would 
determine  the  direction  which  the  project  should  take  from  that  point  on.  The  two  alternatives  would  be  to 
focus  on  the  numerical  solution  of  the  eigenvalue  problem  or  pose  the  problem  of  stability  in  the  open  cavity 
in  the  broader  context  of  solutions  of  (linear)  initial  value  problems.  The  latter  option  would  permit 
addressing  the  issue  of  transient  growth  at  short  times  and  that  of  eigenmodes  at  the  limit  of  large  times  of 
the  same  analysis.  However,  the  cost  of  performing  a  transient-growth  study,  in  addition  to  the  essential 
change  of  focus  in  the  Grant  objectives,  would  require  strong  justification.  It  turned  out  that  linear  theory  in 
the  form  of  BiGlobal  analyses  delivered  unstable  eigenvalues  in  the  model  three-dimensional  problem  and,  as 
such,  the  original  plan  was  pursued. 


Figure  10  illustrates  the  configuration  chosen  in  order  to  address  the  question  of  three-dimensional 
instability  and  relevance  of  BiGlobal  theory  in  the  model  problem.  Motion  of  a  lid  at  an  arbitrary  angle  (f>  sets 
up  a  three-component  basic  flow  inside  the  cavity,  the  limiting  cases  of  which  are  the  classic  lid-driven  cavity 
(at  4>  =  0,  where  BiGlobal  theory  predicts  excellent  agreement  with  experiment;  Theofilis  2000),  and  the  wall- 
bounded  Couette  flow  at  (|)  =  n  /  2  which  is  (erroneously  predicted  by  BiGlobal  analysis  to  be)  stable.  The 
introduction  of  a  third  velocity  component  by  the  lid-motion  along  the  homogeneous  direction  results  in  an 
order-of-magnitude  increase  of  the  computing  effort  in  order  to  solve  the  BiGlobal  EVP;  however,  reliable 
solutions  to  this  problem  can  be  obtained  using  a  purpose-built  extension  of  a  well-validated  incompressible 
BiGlobal  EVP  solver  (Theofilis,  Fedorov,  Obrist  and  Dallmann  2003).  Results  obtained  in  the  3D  cavity 
model  at  all  intermediate  (|) —values,  shown  in  Figure  11,  demonstrate  that  the  known  modes  of  the  lid- 
driven  cavity  S,  Tl,  T2  and  T3  persist  in  the  case  of  three-dimensional  flow.  Three-dimensionality  was 
found  to  be  essential  in  introducing  linear  amplification;  at  the  low  Reynolds  number  examined,  at  which 
steady  two-dimensional  states  were  predicted  by  2d  DNS,  instability  appears  in  the  flow  on  account  of 
P  =  2ti/Lz  ^  0,  Lz  being  the  finite  spanwise  extent  of  the  cavity.  The  levels  of  (viscous)  amplification  (  cy  >  0  ) 
are  analogous  to  those  of  the  lid-driven  cavity,  which  in  turn  agree  very  well  with  experiment.  It  was  further 
discovered  that  the  range  of  amplified  frequencies,  denoted  by  00r  in  Figure  11,  is  a  function  of  the  angle 
considered,  being  widest  at  the  cavity-related  limit.  It  thus  appears  that  the  wall-bounded  Couette  flow  is  a 
singular  exception,  due  to  the  absence  of  velocity  components  on  the  plane  normal  to  the  wavenumber 
vector.  This  result  builds  confidence  in  employing  BiGlobal  analysis  to  the  three-dimensional  compressible 
open  cavity  linear  instability  problem. 
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A  major  publication  (Theofilis,  Duck  and  Owen  2004)  has  thus  been  generated  during  the  first  year  of 
the  Grant.  Material  has  been  composed  of  (a)  BiGlobal  instability  results  in  rectangular  duct  and  wall- 
bounded  Couette  flow,  performed  prior  to  this  grant  in  collaboration  with  Prof.  P.  W.  Duck  and  his  student 
Mr.  J.  Owen  of  Manchester  University,  UK,  (b)  EOARD-supported  research  of  Theofilis  in  the  lid-driven 
cavity  (contract  F61775-99-WE090)  and  (c)  a  subset  of  the  results  obtained  by  nu  Modelling  s  the  first  year  of 
the  present  grant.  Based  on  the  Referees’  report  it  is  clear  that  the  inclusion  of  the  model-cavity  results  in  (c) 
has  been  instrumental  in  the  acceptance  of  the  paper  for  publication. 


Figure  10:  A  model  three-dimensional  flow  in  cavities  (Theofilis.  Duck  &  Owen  J.  Fluid  Mech.  505.  249-286.  2004). 


Figure  11:  Dependence  of  amplification  rates  on  wavenumber  in  (a)  and  frequency  in  (b)  at  three  angles  (|)  =  7t  /  8,  (|)  =  71  /  4 

and  4>  —  3tc  /  8  and  two  Reynolds  numbers  Re  =900  (dotted  line)  and  1000  (solid  line)  (Theofilis.  Duck  &  Owen  J.  Fluid 

Mech.  505.  249-286.  2004i. 


M  +- W 
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5.  Direct  solution  of  the  compressible  BiGlobal  eigenvalue  problem 

Motivated  by  the  need  to  further  understand  whether  and  how  three-dimensionality  changes  the 
transition  boundaries  between  shear  layer  and  wake  mode,  a  new  compressible  BiGlobal  EVP  solver  has 
subsequently  been  designed,  hard-coded  and  validated  at  nu  modelling.  Design  characteristics  of  the  code  have 
been  a)  high-accuracy,  b)  flexibility  in  terms  of  both  the  complexity  of  the  geometry  addressed  and  general 
constitutive  laws  admitted  by  the  code  (equation  of  state,  viscosity  dependence  on  temperature,  the  latter 
permitting  us  to  address  flow  instability  up  to  the  hypersonic  regime)  and,  c)  completeness  and  efficiency  in 
the  computation  of  the  eigenspectrum.  In  the  first  implementation  of  the  code,  steps  taken  towards  meeting 
the  aforementioned  characteristics  have  been  the  use  of  a)  single-domain  spectral  collocation  methods  on  a 
two-dimensional  staggered  Chebyshev  Gauss-Lobatto/ Gauss  grid,  later  to  be  extended  by  spectral- 
multidomain  or  spectral-element  techniques  (Theofilis  EOARD  Final  Report  F61 775-01 -WE049),  while  b)  a 
perfect-gas  law  and  c)  the  QZ  algorithm  for  the  recovery  of  the  entire  eigenspectrum  completed  the  initial 
implementation. 

An  essential  aspect  of  every  new  development  is  careful  validation  against  well-established  results.  The 
compressible  BiGlobal  EVP  as  defined  here  describes  small-amplitude  perturbations  superimposed  upon  an 
arbitrary  two-dimensional  basic  state.  In  the  limiting  case  of  such  a  basic  state  being  a  one-dimensional  profile 
(e.g.  those  pertaining  to  the  boundary-layer  or  shear-layer,  both  of  which  are  directly  relevant  to  the  cavity 
problem  at  specific  locations  in  the  cavity)  the  BiGlobal  two-dimensional  EVP  reduces  to  the  well-known  one- 
dimensional  EVP.  The  latter  EVP  has  been  utilized  to  perform  exhaustive  investigations  into  the  stability  of 
compressible  flat-plate  and  axisymmetric  boundary-layer  and  shear-layer  flows  over  the  span  of  several 
decades  (Mack  1984).  The  newly-developed,  more  general,  code  has  been  validated  thoroughly  against  such 
results,  since  the  latter  form  the  cornerstone  of  practically  all  current  analysis  and  control  methodologies  of 
flow  in  the  open  cavity  that  are  based  on  flow  instability. 

Results  are  presented,  beginning  with  cases  that  serve  to  validate  the  numerical  method  and  illustrate 
our  approach.  These  include  an  incompressible  square  duct,  the  flat-plate  boundary  layer  from 
incompressible  speed  to  Mach  number,  Ma  —  6  flow,  and  an  isolated  parallel  two-dimensional  shear  layer. 
Excellent  agreement  with  the  results  of  Mack  (1984)  has  been  obtained;  as  an  example,  components 
(amplitude  functions  of  disturbance  velocity  and  disturbance  temperature)  of  the  (one-dimensional) 
eigenvector  pertaining  to  three-dimensional  instability  of  a  Ma—  6  flat-plate  boundary  layer  are  shown  in 
Figure  12.  The  complex  structure  of  these  perturbations  and  the  characteristic  of  compressible  flow  critical 
layer  (i.e.  the  location  where  the  phase  velocity  of  the  instability  wave  equals  that  of  the  basic  state,  located  at 
a  wall-normal  distance  y  ~  20  at  these  parameters)  can  clearly  be  seen  in  these  results.  Results  are  presented 
on  the  collocation  grid  used,  in  order  for  the  resolution  requirements  in  the  cavity  problem  to  be  appreciated: 
the  grids  shown  are  the  one-dimensional  restrictions  of  a  two-dimensional  BiGlobal  instability  analysis  grid 
that  would  be  necessary  in  order  to  recover  these  instabilities  at  the  same  level  of  accuracy.  The  amplification 
rates  of  linearly  unstable  compressible  flow  over  a  flat  plate,  accurate  to  all  significant  digits  presented,  are 
shown  in  Table  5;  those  pertinent  to  compressible  shear-layer  flow,  alongside  comparisons  against  earlier 
works,  are  shown  in  Table  6. 
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Figure  12:  Streamwise  and  spanwise  disturbance  velocity  components  and  temperature  perturbations  on  a  flat-plate 
boundary  layer  at  Ma  =  6  (Theofilis  and  Colonius  AIAA  2004-2544). 


In  Table  5  Ma  denotes  Mach  number,  8*  is  the  displacement  thickness  of  the  similarity  solution  of  the 
compressible  flat-plate  boundary  layer  solution  in  multiples  of  a  unit  length,  Re  §*  is  the  associated  Reynolds 
number,  a,  [3  are  wavenumber  parameters  of  three-dimensional  Tollmien-Schlichting  instability  waves,  \\f  the 
associated  angle  (in  radians)  and  (oor ,  00^  the  computed  eigenvalue.  The  same  symbols  are  used  to  describe 
two-dimensional  (P=0)  shear-layer  instability  in  Table  6.  Note  that  the  latter  inviscid  instability  results  have 
been  obtained  by  setting  a  large  Reynolds  number  value,  Re  =  104  ,  in  our  code. 


Table  5:  Validation  of  the  compressible  BiGlobal  EVP  solver  on  the  compressible  flat  plate  boundary  layer  instability 


Ma 

Re  §* 

0 

a 

\|/  =  atan(p/a) 

(CDr  ,  03;) 

10-3 

800 

1.720788 

0.12 

30 

( 0.040583234,  0.001728706) 

1 

1000 

2.004500 

0.10 

45 

(0.038017729,  0.001350618) 

6 

1000 

17.696682 

0.08 

45 

( 0.072487839,  0.001052149) 

Table  6:  Validation  of  the  compressible  BiGlobal  EVP  solver  on  the  compressible  shear  layer  instability 
Ma  a  Blumen  Macaraeg  et  al.  Present  results 


0.0 

0.445 

0.190 

0.18954 

0.189539 

0.2 

0.426 

0.181 

0.18112 

0.181156 

0.4 

0.409 

0.158 

0.15760 

0.157581 

0.6 

0.370 

0.122 

0.12180 

0.121763 

0.8 

0.279 

0.078 

0.07760 

0.077607 

Subsequent  work  at  nu  modelling  has  focused  on  two  aspects  of  the  newly-developed  compressible 
BiGlobal  eigenvalue  problem  (EVP)  solver.  First,  validations  have  continued  using  analytically-known  two- 
dimensional  basic  states.  Results  obtained  independently  by  the  established  incompressible  BiGlobal  EVP 
solver  of  Theofilis  (Theofilis  et  al.  2003;  2004)  have  been  reproduced  by  the  new  compressible  BiGlobal 
solver,  the  latter  run  at  a  low  Mach  number  Ma  —  O(10"4).  One  such  comparison  is  shown  in  Figure  13, 
where  the  excellent  agreement  between  the  eigenmodes  of  rectangular  duct  flow,  as  delivered  by  the  two 
codes  can  be  seen,  alongside  newly-found  instability  modes  associated  with  the  energy  equation;  the  latter  can 
only  be  obtained  by  the  new  compressible  BiGlobal  EVP  solver. 
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Figure  13:  The  BiGlobal  eigenspectrum  of  incompressible  rectangular  duct  flow  (Theofilis,  Duck  &  Owen  2004  /.  Fluid Mech.  505. 

249-286),  recovered  by  the  developed  compressible  BiGlobal  EVP  solver,  using  two  different  two-dimensional  staggered  spectral 

collocation  grids. 
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Second,  the  first-ever  BiGlobal  instability  results  of  three-dimensional  compressible  open  cavity  flow  have 
been  obtained  by  nu  modelling.  The  steady-state  compressible  basic  flow  solution  in  a  L  /  D  =  2  cavity  has 
been  produced  by  2-D  DNS  and  delivered  to  nu  modelling  by  the  Caltech  team.  The  conditions  Re  =  1500  and 
Ma  =  0.325  were  chosen  to  ensure  steadiness  of  the  basic  flow.  In  the  BiGlobal  EVP  solution  the  streamwise 
and  normal  spatial  directions  defining  the  open  cavity  have  been  resolved  in  a  coupled  manner,  while  the 
third,  spanwise  spatial  direction  has  been  treated  as  periodic.  The  extent  of  the  spanwise  domain,  L  }  has  been 
determined  by  a  wavenumber  parameter  fi=2n/Lz  and  chosen  to  be  unity.  First  results  obtained  were  focused 
on  the  flowfield  inside  the  cavity;  at  the  parameters  chosen  only  stable  modes  have  been  recovered.  The 
amplitude  functions  of  the  dominant  streamwise  disturbance  velocity  components  of  the  two  least-damped 
modes  are  shown  in  Figure  14,  where  the  complexity  of  these  structures,  as  well  as  their  connection  with 
instability  mechanisms  outside  the  cavity,  can  clearly  be  seen. 

Figure  14:  The  two  least-damped  BiGlobal  eigenmodes  of  compressible  aspect-ratio-2  open  cavity  flow;  streamwise  disturbance 

velocity  components  of  the  traveling  mode  (left)  and  the  stationary  mode  (right)  (Theofilis  and  Colonius  AIAA  2004-2544) 
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Subsequently,  work  focused  on  two  aspects  of  the  single-domain  compressible  BiGlobal  EVP  solver, 
namely  improving  its  efficiency  and  flexibility.  The  first  improvement  was  accomplished  by  implementing  an 
(implicitly  restarted)  Arnoldi  method  (IRAM)  for  the  solution  of  the  large  partial  derivative  eigenvalue 
problem.  This  algorithm  recovers  the  most  significant  (most  amplified/least  damped)  part  of  the 
eigenspectrum,  as  opposed  to  all  members  of  the  eigenspectrum  that  the  QZ  algorithm  delivers,  as  seen  in 
Figure  15.  The  advantage  of  the  Arnoldi  over  the  QZ  approach  (the  latter  used  by  Theofilis  and  Colonius 
2004)  is  that  the  leading  eigenmodes  are  recovered  using  O (N)  in  the  former,  as  opposed  to  O (N3)  operations 
in  the  latter  algorithm.  N  is  the  leading  dimension  of  the  BiGlobal  eigenvalue  problem  matrix,  which,  in  turn, 
is  the  product  of  nodes  resolving  the  two  spatial  directions  times  the  number  of  equations  solved  (five,  in  the 
compressible  regime).  The  resolution  requirements  are  a  function  of  the  Reynolds  number  addressed;  in 
typical  applications  N  can  be  of  the  order  of  104,  which  renders  the  calculation  speed  advantage  the  Arnoldi 
algorithm  over  the  QZ  approach  decisive  for  the  success  of  BiGlobal  instability  analysis,  in  a  twofold  manner. 

First,  use  of  the  Arnoldi  algorithm  permits  fast  recovery  of  high-resolution  results  at  moderate 
Reynolds  numbers.  Alternatively,  the  Arnoldi  algorithm  delivers  results  of  moderate  (but  adequate)  resolution 
substantially  faster  than  the  QZ  algorithm,  which,  in  turn,  facilitates  performance  of  fast  parametric  studies. 
Second,  use  of  the  Arnoldi  algorithm  requires  storage/manipulation  of  a  single  (large)  array  while  the  QZ 
algorithm  needs  four  such  arrays  to  be  stored.  The  memory  freed  by  the  Arnoldi  algorithm  may  be  devoted  in 
order  to  address  instability  problems  at  high  Reynolds  number  values  inaccessible  to  the  QZ  algorithm.  Using 
the  basic  states  provided  by  Caltech  the  Arnoldi  algorithm  was  employed  to  extend  the  computations  of 
Theofilis  and  Colonius  (2004);  one  result  of  this  effort  is  shown  in  Figure  15.  The  associated  amplitude 
functions  of  the  leading  eigenmodes  are  shown  in  Figure  16  (using  the  QZ  algorithm)  and  in  Figure  17  (using 
the  IRAM  approach).  Both  results  are  obtained  using  the  same  level  of  computing  effort,  explicit  matrix 
formulation  and  storage,  as  well  as  serial  computing.  The  resolution  improvement  that  the  Arnoldi  algorithm 
offers  may  be  appreciated  by  comparing  the  results  of  these  figures. 


Figure  15:  The  BiGlobal  eigenspectrum  of  compressible  two-dimensional  open  cavity  flow  obtained  using  the  QZ  algorithm  (+) 
by  Theofilis  and  Colonius  (2004).  Superimposed  are  shown  identical  values  of  the  most  significant  part  of  the  eigenspectrum.  as 

obtained  using  the  Arnoldi  algorithm  (o). 
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Figure  16:  Amplitude  functions  of  the  disturbance  velocity  components  of  the  first  traveling  (Qr  ^  0)  and  the  first  stationary 

(Qr  =0)  three-dimensional  BiGlobal  eigenmode  at  aspect  ratio  U  /D  =  2  compressible  open  cavity  flow  at  Re  =  1500  and 

Ma  =  0.325,  corresponding  to  a  narrow  spanwise  periodicity  domain  L?  =  L  (]3=27i/Ly  ~  6.28V 


Figure  17:  Amplitude  functions  of  the  disturbance  velocity  components  of  the  first  traveling  (Qr^  0)  and  the  first  stationary 

(Qr  =0)  three-dimensional  BiGlobal  eigenmode  at  aspect  ratio  U/D  =  2  compressible  open  cavity  flow  at  Re  =  1500  and 

Ma  =  0.325,  corresponding  to  a  wide  spanwise  periodicity  domain  L7  =  2,  (fi=2n/L7  ~  3.14V 
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Regarding  flexibility  of  the  BiGlobal  EVP  solver,  results  such  as  those  shown  in  Figures  16  and  17  were 
obtained  by  considering  the  domain  inside  the  cavity  in  isolation  from  the  flowfield  outside  the  cavity  and 
imposing  homogeneous  Neumann  boundary  conditions  in  the  open  end  of  the  domain.  Such  results  were 
confirmed  by  the  independent  Caltech  effort;  this  team  has  extended  their  DNS  code  in  order  to  solve  the 
linearized  Navier-Stokes  equations  and  obtained  leading  eigenmode  results  in  close  qualitative  agreement  with 
those  obtained  by  the  direct  solution  of  the  BiGlobal  eigenvalue  problem;  direct  quantitative  comparisons 
were  not  possible  since  the  computational-resources-intensive  nature  of  the  linearized  Navier-Stokes 
methodology  only  permits  a  modest  resolution  of  the  structures  delivered  by  the  BiGlobal  analysis  inside  the 
cavity.  Conversely,  the  BiGlobal  methodology  at  the  time  of  those  comparisons  was  confined  to  the  solution 
of  the  eigenvalue  problems  in  single-domains,  and  resolved  the  domain  inside  the  cavity,  neglecting  parts  of 
the  eigenmode  potentially  extending  outside  the  cavity,  that  were  threated  by  the  linearized  DNS  approach. 
This  provided  strong  motivation  for  the  developments  that  followed  in  the  direction  of  multi-domain 
solutions  of  eigenvalue  problems,  in  order  to  address  BiGlobal  instability  of  both  the  entire  cavity  domain 
and  complex  full-bay  configurations. 


Figure  18:  The  streamwise  disturbance  velocity  component  of  the  least-damped  BiGlobal  eigenmode  of  compressible 

aspect-ratio-2  open  cavity  flow,  delivered  by  the  Arnoldi  algorithm;  in  the  background  the  steady  laminar  basic  state  upon 

which  this  instability  develops  is  shown. 


M=  0.32 


In  summary,  the  direct  solutions  of  the  BiGlobal  eigenvalue  problem  obtained  in  this  part  of  the 
research  have  been  the  first-ever  successful  solutions  of  the  problem  of  linear  instability  in  the  open  cavity, 
without  the  need  to  resort  to  the  strong  assumption  of  parallel  flow,  as  invariably  done  in  previous  efforts. 
As  such,  structures  inaccessible  to  classic  linear  analysis  of  the  Kelvin-Helmholtz  /  Tollmien-Schlichting  class 
have  been  obtained  for  the  first  time,  as  schematically  depicted  in  Figure  18.  With  work  underway  at  Caltech, 
using  either  the  residuals  algorithm  or  the  linearized  Navier-Stokes  approach  in  order  to  identify  stability 
boundaries  in  the  multiparametric  (empty)  open  cavity  problem,  attention  was  diverted  to  the  issue  of  a  full- 
bay  configuration.  As  an  aside,  it  is  worth  mentioning  that  the  new  single-domain  compressible  BiGlobal 
EVP  solver  was  employed  to  study  another  instability  problem  of  long-standing  interest,  that  in  a  swept 
attachment  line  boundary  layer  in  compressible  flow,  extending  earlier  work  by  Theofilis,  Fedorov,  Obrist 
and  Dallmann  (2003).  Using  the  BiGlobal  EVP  code,  Theofilis  was  able  to  verify  (and  modify)  analytical 
predictions  of  Fedorov  (Moscow  Institute  of  Physics  and  Technology,  RU),  while  Collis  (Sandia  Labs,  USA) 
is  preparing  DNS  work  on  the  infinite  swept  cylinder  geometry,  in  order  to  compare  with  results  of  the 
previous  two  approaches.  Results  were  presented  in  the  invitation-only  IUTAM  Symposium  “One  Hundred 
Years  Boundary  Hayer  Research”,  Gottingen,  Germany,  Aug.  12-14,  2004  and  a  first  publication  on  this  topic  (of 
ongoing  research)  has  been  generated  and  will  appear  in  the  proceedings  of  this  prestigious  Symposium. 
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6.  New  developments  for  open  and  full-bay  models 

6. 1 .  Spectral  multi-domain  algorithms 

One  of  the  main  conclusions  of  the  scientific  exchanges  with  the  Caltech  group  has  been  the 
identification  of  the  need  for  a  novel  multi-domain  discretization  of  the  partial-differential-equation-based 
BiGlobal  eigenvalue  problem,  a  topic  at  the  heart  of  the  project  owing  to  the  geometry  of  two-  and  three- 
dimensional  cavities.  The  need  to  work  on  this  subject  matter  may  be  well  inferred  from  the  results  of  Figures 
14,  16,  17  and  18,  where  it  is  clear  that  the  structures  unraveled  by  the  single-domain  BiGlobal  EVP  solver 
and  confirmed  by  the  linearized  Navier-Stokes  approach  extend  into  the  flow  domain  outside  the  cavity. 
Developments  were  thus  taken  up  by  the  nu  modelling  team  in  order  to  devise  a  new  multi-domain 
compressible  BiGlobal  EVP  solver.  Given  that  the  group  of  Colonius  has  focused  on  high-order  finite- 
difference  schemes  it  was  decided  by  Theofilis  to  work  with  an  alternative  numerical  methodology,  namely 
spectral  collocation.  In  this  manner,  value  was  added  to  the  project  in  terms  of  the  ability  to  cross-validate 
results  using  independent  numerical  methodologies. 

Initial  work  on  the  subject  of  multi-domain  discretization  of  partial  differential  equations  focused  on 
four  aspects:  solution  of  linear  systems  and  eigenvalue  problems  in  one  and  two  spatial  dimensions, 
respectively.  Model  elliptic  problems,  such  as  a  one-dimensional  linear  equation  having  as  solution  the 
quadratic  function  and  the  two-dimensional  linear  equation 

Au  =  10*sin(8*x*(y-l))  in  Q  =  [-1,  1  ]x[-l,  1  ] 

u  =  0.  in  dCl 

were  solved.  Results  of  spectral  multi-domain  solutions  of  the  ID  linear  equation,  defined  on  x  E  [-2,  2],  and 
of  the  2D  problem  defined  above  are  shown  in  Figure  19.  Excellent  performance  of  the  multi-domain 
scheme  may  be  seen  in  the  solution  of  both  the  1-  and  the  2-D  problem. 

Subsequent  work  built  upon  the  success  of  the  previous  efforts  in  solving  elliptic  partial  differential 
equations  and  extended  this  technology  to  tackling  eigenvalue  problems.  Two  classes  of  EVPs  were  solved, 
the  Orr-Sommerfeld  equation,  which  is  the  1-D  limit  of  the  incompressible  BiGlobal  EVP,  as  well  as  partial- 
differential-(Poisson-)-equation-based  scalar  two-dimensional  eigenvalue  problems.  Both  are  essential 
stepping  stones  in  order  to  be  able  to  solve  the  BiGlobal  EVP  in  the  open  cavity.  A  key  validation  step 
included  well-studied  limiting  cases  based  on  the  equations  of  fluid  flow  motion,  the  best  known  of  which  is 
the  Orr-Sommerfeld  equation.  The  latter  equation  has  been  solved  using  two  one-dimensional  domains 
connected  by  imposing  C1  continuity  at  y=0,  and  resolved  using  32  collocation  points  per  subdomain.  The 
parameters  shown  by  Canuto,  Hussaini,  Quarteroni  and  Zang  (1987),  Re=7500,  a=l,  have  been  chosen  for 
comparison.  Excellent  agreement  between  the  corresponding  eigenvalue  and  the  published  results  has  been 
obtained,  while  the  dominant  real  part  of  the  eigenfunction  is  presented  in  graphical  form  in  Figure  20. 

The  excellent  agreement  between  the  literature  results  and  those  delivered  by  our  novel  spectral  multi- 
domain  technique  suggests  that  the  latter  is  well-suited  to  study  the  open  cavity  flow  BiGlobal  instability 
problem.  However,  prior  to  addressing  the  Navier-Stokes  based  EVP,  the  last  validations  performed  have 
concerned  eigenvalue  problems  involving  the  Poisson  equation,  which  is  an  essential  part  of  the  linearized 
operator  describing  BiGlobal  instability  analysis  (Theofilis  2003). 


Figure  19:  Upper:  a  quadratic  one-dimensional  function  obtained  on  four  domains  comprising  8.8.16  and  32  collocation 


34 


FINAL  REPORT 

Grant  FA8655-03-1-3059  (Theofilis)  -  “A  Unified  View  of  Global  Instability  of  Compressible  Flow  over  Open  Cavities” 


points.  Lower:  numerical  solution  of  a  two-dimensional  elliptic  problem  on  a  single  2D  domain  comprising  25x25  points 

and  on  four  2D  domains,  each  comprising  12x12  points. 


x 


Figure  20:  Multi-domain  solution  of  the  Orr-Sommerfeld  eigenvalue  problem  in  plane  channel  flow.  Shown  is  the  real  part  of  the 

amplitude  function  of  the  wall-normal  perturbation  at  Re=7500.  (1=  1 :  domains  connected  at  v=0. 


The  model  EVP  based 
Trefethen  (2000): 


on  the  Poisson  equation  that  has  been  addressed  has  been  suggested  by 


-  A u  +  /(x, y) u  -  X u  e/iQ  =  (-l,l) 

u  =  0  en  SQ 
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Taking  /  =  0  the  eigenfunctions  may  be  obtained  analytically  by  separation  of  variables  and  have  the 

form: 


sin(/fl.(A-  +  l))sin(/clO;  +  l)) 


7Z 


where  kx  and  ky  are  integer  multiples  of  —  .  Consequently  the  eigenvalues  are  also  analytically  known: 


—rii2  +  j2\  i,j=l,2,3, 

4 


making  this  problem  a  good  candidate  to  study  the  properties  of  a  novel  numerical  scheme.  Discretization  of 
the  model  EVP  using  spectral  collocation  and  the  multidomain  technique  results  in  the  EVP  of  the  form 

(LHS)  u  =  A  (RHS)  u 


where  X  is  the  sought  eigenvalue,  LHS  =  (-A  +  / (x,  y ))  and,  here,  RHS  is  the  identity  matrix. 

First,  single-domain  computations  delivered  the  following  results  for  the  four  most  significant 

eigenvalues: 

Table  7:  Single-domain  solution  of  the  model  two-dimensional  EVP  — Au  +/(xvJ  u  =  X  u  for/=0. 


Analytical  result 

Present  numerical  solution 

1 

2-  n2  /  4=  4.93480220054467 

4.934802200544 

2 

5-  n1 14  =12.33700550136169 

12.337005501361 

3 

5-  n1 14  =12.33700550136169 

12.337005501363 

4 

8-  n 2  /4  =19.739208802178717 

19.739208802178 

Taking  f(x,y)  =  exp(20(T_^_l))  only  single-domain  numerical  solutions  of  this  problem  are  known, 
one  of  which  is  presented  by  Trefethen  (2000)  in  terms  of  the  eigenvalues: 

Table  8:  Single-domain  solution  of  the  model  two-dimensional  EVP  — Au  +  f(xLv)  u  =  X  u  for/=  exp  [  20  (y-x-1)  ]. 


Trefethen  (2000) 

Present  numerical  solution 

1 

2.116423652153-  n1  / 4 

5.221491363672 

2 

5.023585398303-  n1 1 4 

12.394945792404 

3 

5.548908101834-  n1  / 4 

13.689074320561 

4 

8.642804449790-  n1 ! 4 

21.323357787175 

The  corresponding  eigenfunctions  are  shown  in  Figure  21. 
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Figure  21;  The  first  four  eigenmodes  of  the  2-D  Poisson  EVP  corresponding  to/fxi/j— exb(y-x-1)*  solved  using  single-domain 

spectral  collocation 


The  same  EVP  is  then  solved  using  the  novel  spectral  multi-domain  technique.  Two  and  four 
subdomains  have  been  considered,  each  being  resolved  by  21 2  collocation  nodes.  C1  continuity  of  the  solution 
across  subdomain  boundaries  has  been  imposed.  Results  have  been  obtained  for  two  cases  of  the  forcing 
function  / (  ),  one  corresponding  to  analytically-  and  the  other  to  only-numerically- known  results;  these 

are  presented  in  the  following  tables: 


Table  9:  Multi-domain  solution  of  the  model  two-dimensional  EVP  — Au  +  /fryd  u  =  X  u  for  f  =  0 


Analytical  value 

Present  solution  - 
2  domains 

Present  solution  - 
4  domains 

1 

4.93480220054467 

4.934802200541 

4.934802200538 

2 

12.33700550136169 

12.337005501359 

12.337005501895 

3 

12.33700550136169 

12.337005501362 

12.337005501900 

4 

19.73920880217871 

19.739208802180 

19.739208803263 

Table  10:  Multi-domain  solution  of  the  model  two-dimensional  EVP  — Au  +/(Vd  u  —  X  u  for  f  —  exp  [  20  (y-x-1)  ]. 


Trefethen  (2000) 
single-domain 

Present  solution 
single- domain 

Present  solution  - 
2  domains 

Present  solution  - 
4  domains 

1 

2.116423652153-  n2 1 4 

5.221491363672 

5.221445153111 

5.221849381051 

2 

5.023585398303-  tt2  1 4 

12.394945792404 

12.394904449197 

12.395128398195 

3 

5.548908101834-  n1  1 4 

13.689074320561 

13.688868122805 

13.690499067793 

4 

8.642804449790 ■  n2  ! 4 

21.323357787175 

21.323199763726 

21.324400372510 

while  the  corresponding  eigenfunctions,  obtained  in  two  domains  joined  at  x=0,  are  shown  in  Figure  22. 
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Figure  22:  The  first  four  eigenmodes  of  the  2-D  Poisson  EVP  solved  using/6cyy] -expfy-x-  /  J.  solved  using  the  novel  multi-domain 

spectral  collocation  technique. 
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The  excellent  agreement  between  analytical  and  numerical  results  has  built  additional  confidence  in  the 
newly-developed  methods  in  order  to  apply  them  to  the  BiGlobal  eigenvalue  problem.  For  the  sake  of 
simplicity  of  the  presentation  of  the  methodology  followed,  in  Figure  23  flow  in  a  square  domain  in  the 
incompressible  limit  is  considered.  The  square  is  subdivided  in  four  subdomains,  denoted  by  the  colors  red, 
green,  dark-  and  light  blue.  The  disturbance  equations  are  written  in  each  of  the  domains,  as  schematically 
shown.  As  a  matter  of  fact,  an  essential  requirement  of  the  spectral  multidomain  algorithm  developed,  not 
shown  in  Figure  23,  has  been  its  non-conforming  nature,  namely  its  ability  to  discretize  neighboring  domains 
with  different  resolutions.  This  requirement  was  imposed  on  the  algorithm  with  a  view  to  be  able  to  resolve 
subdomains  away  from  the  solid  surface  of  the  cavity  by  relatively  coarse  meshes,  while  selectively  reserving 
high  resolution  for  boundary  layers  and  downstream  corners  in  the  cavity.  This  requirement  added  to  the 
complexity  of  the  overall  non-conforming  spectral  multidomain  algorithm,  especially  with  respect  to  the 
imposition  of  distinct  types  of  boundary  conditions.  Figure  24  provides  a  graphical  representation  of  the 
situation  encountered  at  the  interface  of  non-conforming  elements,  while  Figure  25  extends  and  completes 
Figure  23,  in  that  the  means  of  imposition  of  the  two  types  of  Dirichlet  or  Neumann  boundary  conditions,  as 
well  as  that  of  the  compatibility  conditions  at  the  interface  are  shown. 

Figure  23:  Extension  of  the  spectral  multi-domain  technique  to  the  BiGlobal  EVP.  Conforming-domain  discretization. 
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Figure  24:  Typical  interface  of  two  subdomains  in  the  non-conforming  spectral  multidomain  technique  developed. 


Figure  25:  Details  of  the  imposition  of  Dirichlet-  and  Neumann- type  of  boundary  conditions  and  of  compatibility  conditions  at 

the  interface  of  neighboring  subdomains  in  the  solution  of  the  BiGlobal  EVP  by  the  non-conforming  spectral  multidomain 

technique. 
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6.2.  Yinite-element  algorithms 

In  view  of  the  ultimate  aim  of  this  avenue  of  research,  instability  analysis  of  complex  full-bay  models  of 
open  cavity  flows,  the  necessity  arises  to  compute  steady  basic  states  around  such  configurations.  One  option 
would  have  been  to  extend  and  adapt  the  Caltech  high-order  accurate  DNS  code  to  provide  such  states  on 
conformable  regular  (rectangular)  subdomains.  Instead,  we  opted  for  finite-element  solutions  on  unstructured 
meshes,  an  approach  which  permits  resolution  of  both  structured  rectangular  domains,  as  well  as  more 
complex  irregular  geometries  of  cavity  bays,  and  ultimate  consideration  of  arbitrarily-shaped  (including 
curved)  stores.  During  the  current  validation  stages,  assumptions  were  relaxed  in  a  systematic  manner, 
departing  from  the  well-understood  and  documented  rectangular  cavity  configuration.  The  cavity  geometry 
and  unit  Reynolds  number  have  been  defined,  scaling  lengths  with  the  cavity  depth,  D.  An  analytic  model  of 
the  (Blasius)  flat  plate  boundary  layer  has  been  developed  and  imposed  at  the  inflow  boundary.  In  order  to 
recover  the  steady  basic  states  to  be  analyzed,  points  in  the  multi-parametric  space  (Reynolds  number, 
thickness  of  the  boundary  layer  at  the  cavity  lip,  length-over-depth  of  the  cavity)  have  been  identified,  at 
which  such  steady  states  exist.  Basic  flows  have  been  computed  on  unstructured  grids,  typically  comprising 
O(105)  triangular  elements,  distributed  in  a  manner  that  refines  critical  calculation  areas  such  as  the  boundary 
layers  and  the  cavity  corners;  a  typical  mesh  at  Re  =1600  is  shown  in  Figure  26. 

Figure  26:  Typical  unstructured  mesh  used  for  the  finite  element  solution  of  the  basic  flow  problem,  here  at  Re=1600. 


Regarding  the  spectral  multi-domain  BiGlobal  eigenvalue  problem  solution  algoritm  subsequently  used, 
the  basic  flow  obtained  by  the  finite-element  method  was  interpolated  from  a  subset  of  the  unstructured 
mesh  on  each  subdomain  of  the  spectral  multi-domain  grid,  in  a  manner  that  will  be  described  shortly.  In 
order  to  resolve  well  areas  where  the  sought  eigenvectors  are  expected  to  attain  their  maxima,  two  strategies 
have  been  investigated.  First,  the  geometric  features  of  the  open  cavity  were  taken  into  account,  dividing 
space  into  four  major  rectangular  subdomains,  one  being  the  cavity  itself,  and  the  other  three  being  defined 
upstream-,  above-  and  downstream  of  the  cavity.  Each  of  the  four  subdomains  was  stored  on  a  different 
processor  of  the  computing  cluster  and  was  discretized  independently,  adding  points  to  the  appropriate 
subdomain  (typically  that  inside  the  cavity)  until  convergence  was  achieved.  This  strategy  was  soon 
abandoned  since  the  resolution  requirements  in  specific  subdomains  are  such  that  the  available  memory  of 
the  computing  node  at  which  those  subdomains  resided  quickly  became  inadequate,  while  the  availability  of 
the  cluster  computing  nodes  was  not  fully  exploited. 
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A  second  strategy  was  then  followed,  in  which  the  same  four  subdomains  have  been  further  subdivided 
and  the  resulting  subdomains  have  been  distributed,  one  per  processor.  Since  the  resulting  subdomains  are 
substantially  smaller  in  size  in  the  latter  strategy,  it  is  straightforward  to  refine  only  those  parts  of  space  where 
the  eigenmode  resides;  in  turn  full  advantage  is  taken  of  the  existence  of  several  processors  of  the  computing 
cluster.  It  has  to  be  noted  that  the  latter  approach  relies  precisely  on  the  availability  of  the  computing  cluster 
described  in  Section  IV.l  and  further  justifies  the  initially  projected  need  for  investment  on  hardware 
appropriate  to  support  the  numerical  work  of  this  Grant.  Parallelism  has  been  achieved  through  performing 
the  residual  matrix  operations  using  iterative  algorithms,  as  implemented  in  standard  (distributed)  libraries. 

The  overall  computational  strategy  for  the  recovery  of  BiGlobal  instability  results  using  the  non- 
conforming  spectral  multidomain  algorithm  is  shown  in  Figure  27.  A  key  difference  with  the  previously  used 
QZ  and  IRAM  approaches  for  the  computation  of  the  eigenspectrum  is  that  in  the  present  distributed 
computing  strategy  the  matrix  discretizing  the  eigenvalue  problem  is  never  formed  explicitly,  although  the 
(parallel)  Arnoldi  algorithm  is  still  used.  By  contrast  to  the  previously  used  approach,  the  matrix  is  formed 
using  the  compressed  column  format  (CCF)  scheme,  as  opposed  to  the  dense  matrix  storage  used  in  the 
BiGlobal  solutions  presented  in  Sections  IV.  4  —  5.  In  the  present  CCF  formulation  only  the  non-zero  entries 
of  the  matrix  are  kept  in  memory  and  are  then  processed  by  the  iterative  solution  software. 


Figure  27:  Overall  solution  algorithm  for  the  BiGlobal  EVP  resulting  from  non- conforming  spectral  multidomain  spatial 

decomposition. 
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With  all  the  algorithms  in  place  and  validated,  efforts  were  then  directed  toward  the  solution  of  the 
incompressible  BiGlobal  EVP  in  a  simple  geometry,  previously  analyzed  in  the  EOARD-funded  work  of 
Theofilis  et  al.  (2004),  namely  stable  flow  in  the  lid-driven  cavity  at  Re  =200,  (3=1.  Both  the  single-domain 
case  and  two  multi-domain  cases  were  solved  with  the  same  newly  developed  code.  The  minimum  nontrivial 
discretization  involves  partition  of  the  square  into  two  sub-domains,  denoted  “red”  and  “green”,  each  of 
which  is  discretized  using  the  same  number  of  collocation  points  as  in  the  single-domain  case,  32  in  each 
spatial  direction.  Excellent  quantitative  agreement  has  been  obtained  between  the  literature  value  and  those 
obtained  using  the  multi-domain  technique,  maximum  differences  in  the  eigenvalues  being  confined  below 
the  fourth  significant  place.  The  amplitude  functions  of  the  leading  member  of  the  eigenspectrum  obtained 
through  the  three  approaches  are  shown  in  Figure  28.  These  results  demonstrate  that  the  first-ever 
solution  of  the  incompressible  BiGlobal  eigenvalue  problem  in  multiple  domains  has  been  obtained 
in  the  framework  of  the  present  Grant. 


Figure  28:  The  spatial  structure  of  the  streamwise-  and  transverse-  disturbance  velocity  and  disturbance  pressure  amplitude 

functions  pertaining  to  the  leading  eigenvector  of  the  lid-driven  cavity  flow  at  Re=200,  P=h  as  delivered  by  single-domain 

and  two  alternative  multi-domain  solutions  of  the  BiGlobal  EVP  (Theofilis.  Duck  &  Owen  TFM  505  2004) 


Single-domain  BiGlobal  analysis 
(Theofilis,  Duck  &  Owen 
JFM  505,  249-286,  2004) 


Multi- domain  BiGlobal  analysis 
2  domains  connected  at  x=0.5 
(present  analysis) 


Multi- domain  BiGlobal  analysis 
2  domains  connected  at  y=0.5 
(present  analysis) 


42 


FINAL  REPORT 

Grant  FA8655-03-1-3059  (Theofilis)  -  “A  Unified  View  of  Global  Instability  of  Compressible  Flow  over  Open  Cavities” 


Next,  attention  has  been  paid  to  the  open  cavity  configuration.  A  steady  laminar  basic  flow  was 
obtained  by  the  finite-element  algorithm  at  Re=1600,  as  shown  in  Figure  29.  The  marked  area  of  the 
flowfield  in  the  neighborhood  of  the  cavity  was  extracted  and  interpolated  on  a  non-conforming  spectral 
multidomain  grid,  as  also  shown  in  this  figure.  Grid  independence  tests  ensured  the  quality  of  the 
interpolations  and  provided  the  variable  coefficients  for  the  subsequent  numerical  solution  of  the  three- 
dimensional  BiGlobal  eigenvalue  problem. 

Figure  29:  Left  column  shows  typical  basic  flows  obtained  by  the  finite-element  method.  The  area  within  the  dotted  box  is 

extracted  and  analyzed.  Right  column:  close-up  of  the  interpolated  basic  flow  velocity  components,  with  the  non-conforming 

spectral  multi-domain  discretization  clearly  visible. 


Original  basic  flow  streamwise  (upper)  and 
transverse  (lower)  velocity  components 


Interpolated  basic  flow  velocity  components 


The  first  two  stationary  modes,  obtained  at  Re  =1600,  and  spanwise  wavenumber  /?  =  1,  are  shown  in 
Figure  30  in  terms  of  their  respective  three  components  of  the  disturbance  velocity  amplitude  functions.  Here 
16  subdomains  have  been  used,  those  outside  the  cavity  having  been  resolved  using  12  spectral  collocation 
points  in  each  spatial  direction,  while  each  of  the  four  subdomains  inside  the  cavity  has  been  resolved  using 
24  spectral  collocation  points  per  direction.  It  is  interesting  to  notice  that  at  these  parameters,  instability 
activity  is  confined  within  the  cavity,  where  all  amplitude  functions  attain  their  respective  maxima,  while  the 
disturbance  flowfield  outside  the  cavity  remains  benign.  These  conditions  are  reminiscent  of  the  instability 
mechanisms  unraveled  by  the  residuals  algorithm  and  presented  in  Figure  8,  both  situations  being  apparently 
related  with  the  wake  mode  instability.  In  such  situations,  a-posteriori  justification  for  the  (computationally 
less  intensive)  approach  of  Theofilis  and  Colonius  (2004)  may  be  provided.  However,  other  types  of 
disturbances  extending  beyond  the  cavity  may  well  be  associated  with  this  configuration,  which  can  only  be 
resolved  by  the  present  general  approach.  We  refrain  from  further  (speculative)  elaboration  on  these  point 
here  but  note  that  linear  amplification  of  such  modes  will  lead  the  steady  laminar  basic  flow  shown  in  figure 
29  into  unsteadiness  (on  account  of  growth  of  the  eigenmode  and/or  its  nonzero  frequency)  and  three- 
dimensionality  (on  account  of  /?  =/  0)  and  it  is  worth  investigating  whether  turbulence  or  nonlinear  limit  cycles 
is  the  result  of  such  amplification.  These  questions,  as  well  as  further  physical  interpretation  of  these  results, 
are  the  topics  of  a  publication  in  preparation. 
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Figure  30:  The  first  three-dimensional  solutions  of  the  direct  BiGlobal  EVP  in  an  open  cavity  configuration  at  Re— 1600,  (3=1 

incorporating  the  entire  cavity  domain  in  the  solution.  Shown  are  the  first  two  stationary  modes  in  the  cavity. 


However,  worthy  of  mention  in  these  results,  and  others  obtained  but  not  shown  here  in  the  same 
parameter  range,  is  the  qualitative  analogy  of  the  spatial  structure  of  the  leading  stationary 
eigenmodes  of  the  open  and  the  lid-driven  cavity  flow  (Theofilis  2000).  For  comparison,  the  spanwise 
disturbance  velocity  component  of  the  lid-driven  cavity  flow  is  shown  in  Figure  31.  It  is  interesting  that 
structures  corresponding  to  essentially  nonparallel  flows  are  present  in  both  cavity  configurations.  None  of 
them  appears  to  be  related  with  the  well-studied  KH  shear-layer  instability,  which  is  not  surprising  in  the  lid- 
driven  but  is  remarkable  in  the  open-cavity  configuration. 

Figure  31:  The  spanwise  disturbance  velocity  component  of  the  lid-driven  cavity  flow  Theofilis  Duck  &  Owen  JFM  505  2004V 
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It  is  a  straightforward  matter  to  analyze  different  open  cavity  configurations  using  the  developed 
methodology.  From  a  physical  point  of  view  it  is  necessary  to  lower  the  Reynolds  number  below  that  used  in 
the  square  cavity  analysis  in  order  to  obtain  a  steady  state,  a  situation  which  is  consistent  with  the  linear 
BiGlobal  instability  scenario  in  the  lid-driven  cavity  (Theofilis  et  al.  2004).  In  Figure  32  results  are  presented 
of  an  analysis  of  L  /  D  =  2  flow  at  Re  =  400  and  P  =  1.  The  leading  travelling  and  stationary  eigenmodes  are 
shown  in  terms  of  the  magnitude  of  the  amplitude  functions  of  the  three  disturbance  velocity  components. 
Interestingly,  the  (spanwise-)  travelling  disturbance  is  also  situated  inside  the  cavity,  while  the  stationary  mode 
extends  in  the  shear-layer  region.  The  three-dimensional  (out-of-plane)  disturbance  velocity  component  of 
both  modes  shows  strong  peaks  in  the  downstream  corner  of  the  cavity.  Systematic  studies  of  the 
multidimensional  parameter  space  are  currently  underway;  however,  the  following  discussion  focuses  on  a 
different  aspect  of  cavity  instability,  using  the  L/D=2  results  of  Figure  32  as  reference. 


Figure  32:  Amplitude  functions  of  the  leading  traveling  and  stationary  mode  in  an  aspect  ratio  2  open  cavity  at  Re— 400.  spanwise 

wavenumber.  (3=1. 
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6.3.  On  open-  v.  full-bay  cavity  instability 

The  key  question  which  led  to  the  development  of  the  non-conforming  spectral  multidomain  algorithm 
has  been  the  investigation  of  the  potential  influence  that  stores  inside  an  open  cavity  may  have  on  the 
instability  characteristics  of  the  cavity.  Answering  this  question  has  been  facilitated  by  the  newly  developed 
multi-domain  algorithm,  which  permits  defining  and  analyzing  efficiently  any  regular  decomposable  domain, 
concretely  in  the  case  of  the  open  cavity,  any  rectangular  decomposable  domain.  The  simplest  model  store 
configuration  possible,  namely  a  rectangular  object  placed  inside  an  open  cavity,  was  considered  as  a 
demonstrator  and  has  been  analyzed  with  respect  to  its  BiGlobal  instability.  The  basic  states  corresponding  to 
the  aforementioned  model  problem  are  shown  in  Figure  33.  While  the  streamwise  velocity  component  is 
hardly  affected  by  the  presence  of  the  object,  strong  modifications  of  the  transverse  basic  flow  velocity 
component  may  be  seen.  The  main  result  of  the  subsequent  instability  analysis  is  that  this  0(1)  modification 
of  the  basic  state  results  in  significant  departure  of  the  instability  characteristics  of  a  cavity  containing  objects 
from  those  of  the  reasonably  well-understood  (empty)  open  cavity  configuration.  A  demonstration  of  this 
fact  may  be  seen  in  Figure  34,  in  terms  of  the  leading  stationary  eigenmode  of  the  configurations  presented  in 
Figure  32,  at  the  same  flow  parameters. 
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Figuqqre  33:  Upper  row:  Streamwise-  (left)  and  transverse  (right)  basic  flow  velocity  components  of  incompressible  flow  over  an 

aspect  ratio  2  (empty)  open  cavity  at  Re=400.  Lower  row:  Same  result  with  a  model  store  placed  inside  the  cavity. 
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Figure  34:  Upper  row:  Streamwise-  (left),  transverse  (middle)  and  spanwise  (right)  disturbance  flow  velocity  components  of  the 

open  cavity  basic  flow  shown  in  Figure  L  at  Re=400.  spanwise  wavenumber,  (3=1.  Lower  row:  Corresponding  result  with  an  object 

placed  inside  the  cavity. 
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The  final  class  of  problem  solved  involved  full  utilization  of  the  potential  of  the  multidomain  algorithm 
in  terms  of  a  model  store  of  higher  complexity  compared  with  that  shown  in  Figures  33  and  34.  This  so- 
called  “T-store”  configuration,  alongside  the  non-conforming  multidomain  grid  used  to  analyzed  it,  is  shown 
in  Figure  35.  Note  that  29  sub-domains  have  been  used  in  order  to  resolve  the  geometry  properly,  devoting  a 
large  number  of  Chebyshev  Gauss-Lobatto  nodes  in  the  neighborhood  of  large  gradients  and  progressively 
relaxing  the  resolution  requirements  away  from  the  boundary  layers  and  the  cavity/ store  walls.  The  flexibility 
of  the  algorithm  developed  is  borne  out  by  the  non-conforming  nature  of  the  sub-domains.  Its  efficiency  may 
be  appreciated  in  noting  that  each  BiGlobal  eigenvalue  computation  of  this  configuration  at  one  set  of 
parameters  requires  approximately  800  seconds  of  (parallel)  computation.  For  the  analysis  the  same 
parameters  as  those  used  in  the  results  presented  in  Figures  32  and  34  have  been  used. 

Figure  35:  Domain  decomposition  and  2D-Chebvshev  Gauss-Lobatto  grids  utilized  for  the  BiGlobal  instability  analysis  of 

the  T- store  configuration 


Figure  36:  The  leading  BiGlobal  eigenmodes  of  the  T-store  configuration  at  the  parameters  shown  in  Figure  34, 
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Again,  the  particular  characteristics  of  the  T-store  basic  flow  (not  shown  here)  have  been  found  to  have 
a  profound  influence  on  the  corresponding  instability  characteristics.  The  leading  two  eigenmodes  of  the 
open  cavity  in  the  presence  of  the  T-store  model  are  shown  in  Figure  36,  in  terms  of  the  amplitude  functions 
of  the  disturbance  velocity  components.  Besides  the  fact  that  frequencies  of  traveling  modes  are  0(1)  apart 
compared  to  those  found  in  the  analyses  summarized  in  Figure  32  and  34,  it  is  noted  that  the  introduction  of 
additional  corners  in  the  geometry  results  in  additional  peaks  appearing  in  the  amplitude  functions  of  the 
components  of  the  stationary  eigenmodes,  pointing  to  regions  that  potentially  would  need  to  be  controlled  in 
the  flow. 
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Summarizing  this  part  of  the  research,  it  becomes  clear  that  instability  characteristics  such  as 
amplification  rates  and  frequencies  of  the  leading  eigenmodes  in  the  three  open  cavity  configurations 
addressed  (empty,  rectangular-object,  T-store)  are  distinct.  While  from  a  theoretical  point  of  view  either  of 
the  three  (or  indeed  any  other  related)  models  could  be  selected  and  analyzed  in  detail,  our  work  has 
demonstrated  for  the  first  time  that,  if  answers  relevant  to  Air  Force  applications  are  sought,  one  should 
focus  on  a  model  of  an  open-cavity-plus-stores  as  close  as  possible  to  the  actual  application.  In  this  respect, 
the  algorithms  developed  during  the  course  of  the  present  Grant  provide  a  flexible  and  efficient  tool  to 
support  such  efforts. 
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V.  SUMMARY  AND  FUTURE  DIRECTIONS 

The  novel  BiGlobal  instability  analysis  approach  was  applied  for  the  first  time  to  incompressible  and 
compressible  flows  over  open  cavities.  Two  distinct  methodologies  have  been  developed,  an  extension  of  the 
residuals  algorithm,  previously  demonstrated  in  the  lid-driven  cavity  problem,  and  direct  solutions  of  the 
BiGlobal  EVP.  The  latter  methodology  itself  has  followed  two  different  paths,  single-domain  computations 
of  the  eigenspectrum  disturbances  inside  the  cavity  in  isolation  from  the  flowfield  outside  the  cavity,  and 
multidomain  computations  of  the  eigenspectrum  corresponding  to  the  entire  two-dimensional  cavity 
flowfield.  Finally,  the  multidomain  algorithm  itself  was  applied  to  both  (empty)  open  and  full-bay  model 
cavity  configurations. 

Consistently,  the  results  obtained  by  all  novel  approaches  presented  herein  have  opened  new  research 
avenues.  The  residuals  algorithm  offered  a  new  unified  perspective  of  known  hydrodynamic  and  aeroacoustic 
instabilities,  in  addition  to  providing  access  to  previously  unknown  (or  not  recognized  as  such)  linear 
disturbances  of  this  flow;  results  of  the  systematic  application  of  the  residuals  algorithm  will  be  presented  by 
the  Caltech  team.  The  direct  solutions  of  the  BiGlobal  eigenvalue  problem  have  been  the  first-ever  successful 
solutions  of  the  problem  of  linear  instability  in  the  open  cavity,  without  the  need  to  resort  to  the  strong 
assumption  of  parallel  or  near-parallel  flow,  as  invariably  done  in  previous  research.  As  such,  structures 
inaccessible  to  classic  linear  analysis  of  the  Kelvin-Helmholtz  /  Tollmien-Schlichting  /  Parabolized  Stability 
Equations  class  have  been  unraveled  by  the  BiGlobal  approach.  In  view  of  the  BiGlobal  instability  results 
obtained  in  the  framework  of  the  present  effort,  it  is  clear  that  future  instability  analyses  of  cavity  flows, 
aiming  at  producing  a  full  picture  of  linear  mechanisms  at  specific  parameter  ranges,  must  consider  the 
problem  of  small-amplitude  disturbance  excitation  from  the  broader,  BiGlobal,  perspective.  In  turn,  efforts 
aiming  at  providing  theoretically-founded  approaches  to  flow  control  by  control  of  flow  instability  must 
answer  the  question  of  how  to  control  unstable  traveling  and  stationary  BiGlobal  eigenmodes,  in  addition  to 
controlling  well-understood  feedback-loop  mechanisms.  Finally,  if  the  relevance  of  cavity  research  to  the 
warfighter  is  to  be  increased,  the  present  research  has  shown  that  a  close  model  of  the  cavity  geometry, 
including  modeling  the  stores,  is  not  only  necessary  but  also  feasible,  thanks  to  the  numerical  methodologies 
developed  herein. 

Regarding  future  directions,  in  the  opinion  of  the  author,  two  areas  are  considered  as  potentially 
extremely  fruitful.  First,  secondary  BiGlobal  theory,  extending  the  tools  presented  herein  to  include  Floquet 
theory,  could  be  used  in  order  to  analyze  time-periodic  cavity  flows  at  Reynolds  numbers  above  the  relatively 
low  ones  at  which  steady  states  exist.  Second,  the  BiGlobal  approach  could  be  combined  with  Reynolds- 
averaged  Navier-Stokes  (RaNS)  predictions,  in  order  to  analyze  (BiGlobal)  instability  of  the  turbulent  flow  at 
flight  Reynolds  and  Mach  numbers.  Both  suggested  avenues  are  feasible:  secondary  BiGlobal  theory  has 
recently  been  demonstrated  by  the  author  in  the  Low  Pressure  Turbine  instability  problem  (see  Final  Report 
of  Grant  F49620-03-1-0295  on  the  subject).  On  the  other  hand,  Crouch  (2005)  has  arrived  at  excellent 
predictions  of  the  wing-buffeting  problem,  analyzing  for  the  first  time  a  RaNS  solution  by  a  compressible 
BiGlobal  instability  approach.  Both  research  avenues  could  directly  incorporate  the  newly-developed  spectral 
multidomain  techniques  in  order  to  analyze  instability  of  full-bay  configurations  at  realistic  flight  conditions 
and  thus  provide  theoretically-founded  answers  to  issues  directly  relevant  to  the  warfighter. 
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VII.  OTHER  ACTIVITIES  SUPPORTED  BY  GRANT 
1.  Collaborations 

On  account  of  knowledge  attained  in  the  framework  of  the  present  Grant,  the  following  new 
collaborations  have  been  started: 

•  Dr.  S.  Collis  (Sandia  Labs)  and  Dr.  Alexander  Fedorov  (Moscow  Institute  of  Physics  and  Technology) 
on  compressible  attachment-line  instability. 

•  Dr.  Emmanuel  Leriche  (Ecole  Polytechnique  Federale  de  Lausanne,  Switzerland)  on  cavity  flows. 

•  Prof.  Rama  Govindarajan  (Jawarhalal  Nehru  Center  for  Advanced  Scientific  Research,  Bangalore, 
India)  on  large-scale  computations  for  the  BiGlobal  eigenvalue  problem. 

•  Dr.  Olaf  Marxen  (KTH,  Stockholm,  Sweden)  on  numerical  methods  for  the  BiGlobal  eigenvalue 
problem  and  DNS. 

Regular  contacts  are  kept  with  all  these  teams. 


2.  Visits/Student  Exchanges 

Theofilis  spent  short  visits  at  Caltech  (of  typically  1  to  2  weeks)  in  2002,  2003  and  2004.  Partially 
supported  by  the  Grant,  Theofilis  visited  Prof.  I.  Gursul  (Univ  of  Bath,  UK)  and  attended  three  invitation- 
only  meetings,  namely  the  IUTAM  Symposium  “One  hundred  years  boundary  layer  research”,  Goettingen, 
Germany,  Aug.  2004,  the  NATO  RTO  AVT-111  meeting,  Prague,  Czech  Republic,  Nov.  2004,  and  the 
IUTAM  Laminar-turbulent  transition  symposium  VI,  Bangalore  India,  Dec.  2004.  He  plans  to  visit  Dr. 
Leriche  (EPFL)  during  two  weeks  in  July  2006. 

De  Vicente  spent  three  months  in  the  summer  of  2004  at  Caltech  in  the  group  of  Colonius  and 
attended  the  VKI  Short  Lecture  Series  “High-order  methods  for  fluid  dynamics”,  Nov.  2005.  He  plans  to 
visit  Dr.  Leriche  (EPFL)  during  the  month  of  July  2006. 


VIII.  PRESENTATIONS  AND  PUBLICATIONS 
1.  Presentations 

Full v-  or  partially  (in  addition  to  WPS  or  other  funds)  supported  by  the  Grant 

Theofilis  delivered  talks  ate  the  following  meetings: 

•  AFOSR  Contractors5  and  Grantees’  Meetings: 
o  Santa  Monica,  CA  (2002), 
o  Destin,  FL  (2003), 
o  Denver,  CO  (2004), 
o  Long  Beach,  CA  (2005) 
o  Atlanta,  GA  (2006)  —  planned. 
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•  Invited  talks: 

o  UC  San  Diego,  CA,  visiting  Prof.  T.  Bewley 
o  U  of  Arizona,  Tucson,  AZ,  visiting  Prof.  I.  Wygnanski. 

•  AIAA  contributions: 

o  Orlando,  FL  (2003), 
o  Portland,  OR  (2004) 
o  San  Francisco,  CA  (2006)  —  planned. 

De  Vicente  delivered  talks  at  the  following  meetings: 

•  Universidad  Politecnica  de  Madrid.  ''Implementation  de  metodos  numericos  en  maquinas  de  paralelo 
para  problemas  fluidodinamicos”.  February  2006,  Madrid,  Spain. 

•  3rd  Symposium  on  Global  Flow  Instability  and  Control.  "Global  instability  of  realistic  open  cavity 
flows”  September  2005  Crete,  Greece. 

2.  Publications 

Fully  supported  by  the  Grant 

Theofilis,  V.  and  Colonius,  T.  2003  An  algorithm  for  the  recovery  of  2-  and  3-D  BiGlobal  instabilities  of 
compressible  flow  over  2-D  open  cavities.  AIAA  Paper  2004—4143. 

Theofilis,  V.  and  Colonius,  T.  2004  Three-dimensional  instabilities  of  compressible  flow  over  open  cavities: 

direct  solution  of  the  BiGlobal  eigenvalue  problem.  AIAA  Paper  2004—2544. 

Theofilis,  V.,  Duck,  P.  W.,  and  Owen,  J.  2004  Viscous  linear  stability  analysis  of  rectangular  duct  and  cavity 
flows  Journal  of  Fluid  Mechanics  505:249—286. 

Theofilis,  V.,  Fedorov,  A.  and  Collis,  S.  S.  2004  Leading-edge  boundary  layer:  Prandtl’s  vision,  recent 
developments  and  future  directions.  IUTAM  Symposium  "One  Hundred  Years  Boundary  Layer 
Research”,  Gottingen,  Aug.  12-14,  2004,  Springer,  ISBN  1-4020-4149-7. 

De  Vicente,  J.,  Valero,  E.,  Gonzalez,  L.  and  Theofilis,  V.  2006  Spectral  Multi-Domain  Methods  for  BiGlobal 
Instability  Analysis  of  Complex  Flows  over  Open  Cavity  Configurations.  36th  AIAA  Fluid  Dynamics 
Conference  and  Exhibit,  San  Francisco,  CA,  June  5-8,  2006. 

De  Vicente,  J.  and  Theofilis,  V.  2006  Non  conforming  spectral  multidomain  methods  for  BiGlobal  instability 
analysis,  (in  preparation). 

Fartially  supported  by  the  Grant 

Theofilis,  V.,  Seifert,  A.,  Joslin,  R.  D.,  and  Collis,  S.  S.  2004  Poster  on  Theoretical  and  Experimental  Active 
Flow  Control.  AIAA  Summer  Meeting,  Portland,  OR,  June  28  —  July  1,  2004. 

Collis,  S.  S.,  Joslin,  R.  D.,  Seifert,  A.  and  Theofilis,  V.  2004  Issues  in  active  flow  control:  theory,  control, 
simulation  and  experiments.  Progress  in  Aerospace  Sciences  40:237—289.  (Top25 article  in  Prog.  Aero. 
Sci.) 

Theofilis,  V.,  Fedorov,  A.  and  Collis,  S.  S.  Leading-edge  boundary-layer:  Prandtl’s  vision,  recent 
developments  and  future  directions.  IUTAM  Symposium  "One  hundred  years  of  boundary  layer 
research”,  Gottingen,  Germany,  Aug  13-17,  2004.  Springer. 

Simens,  M.,  Gonzalez,  L.,  Theofilis,  V.,  and  Gomez-Bianco,  R.  2005  On  fundamental  instability  mechanisms 
of  nominally  2-D  separation  bubbles.  IUTAM  Laminar-Turbulent  Symposium  VI,  Bangalore,  India, 
Dec.  13-17,  2004.  (R.  Narasimha  and  R.  Govindarajan  eds.)  Springer  2005. 


53 


FINAL  REPORT 

Grant  FA8655-03-1-3059  (Theofilis)  -  “A  Unified  View  of  Global  Instability  of  Compressible  Flow  over  Open  Cavities” 


IX.  APPENDIX 

The  equations  governing  compressible  BiGlobal  linear  instability 

The  linear  operators  appearing  in  the  compressible  BiGlobal  eigenvalue  problem  are  defined  next. 
Over-barred  quantities  denote  basic  flow  and  subscripts  denote  partial  differentiation  with  respect  to  the 
variable  shown.  The  equation  of  state  y  M2  p  =  p  T  has  been  used  explicitly,  while  M  (=Ma),  Pr,  and  Re, 
respectively  denote  Mach,  Prandtl  and  Reynolds  numbers. 
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Disturbance  ^-momentum: 

G&  =  ~T~R^P  ~  ^7”  (^)  Tx+p  ^ 

=  ~k~ikP  ®y  ~  G  +  P 

=  “TV  ^  ^  “  TV  ^  +  [“TV  (^)  T*  +  P  u  ]  ^  +  [“TV  (S)  ^  u  ]  ^ 

-h^7  /i  ^  4-  pux  +pvy  +  i  3  pw  +  u  px  -h  v  py 

w  o  -  3-  ( w  & 

-  ^  He  I  tfr  y  y  y 

H-  ^-?yy  3  1  1 3  3~  ~r7  (  d^T2  )  [^"!a!  ^a’  ^4/  ^y]  "T  ^  ^y) 

&zp  ~  3  ;J/fi  +  +  V  Viy) 

IZzd  =  1  p 


r  .  _ _ 1_  { dp_ 

Re  \(£T 

_ 1_  { 

Re  {dT 


54 


FINAL  REPORT 

Grant  FA8655-03-1-3059  (Theofilis)  -  “A  Unified  View  of  Global  Instability  of  Compressible  Flow  over  Open  Cavities” 


Disturbance  Energy: 


li  =  [-2  1  ^  —  ft  (|  Ux  -  |  Vyj  +  7  P  T  j  dx  -2  ^  £  —  ft  (uy  +  Ifc)  dy 

(  —  j  “i!  +  3  )  +  7  P  ^  J 


“  fit1 


—2  1  ^  T  ^  —  P  WX  +  T  px  +  p Tx 
_2  _£ — L-e p,  4.  dx  4-  |^—2  — 

-2  1  A?°  P  ™v  +  f  Py+P  % 


n -1?  -Ml  n 

Re  ** 


=  -2  7  M2  p  «L  0,  -  2  Xii4^  a  WV  ft 


—  ^  - ]£7 -  (If®  t/a;  —  i  -  ft  UJy 

+1  v  '  1  •  P  (“1  +  v*)  +  i  Pi  P'T 


"  Re  Pr  -  w 
^  T  * 


91-2 


'  fie-  Pr  fie-  Pr 
T  fT-1)  Ml 


_ X _ 

fie  Pr 
d-K  1 
f£T 


(#)  tx  dx 


fll-2 


tffA  J.J,  fle,  p,,  -■“  Jjp, 

(#)  (t„+tw)-te?[p?  (^)  (r2  +  **) 


£ep  =  ?j.  ^  +  v  dy  +  i  /?  w  +  7  +  tty  ) 


(^)  {  1  Vl  +  ^  +  I  Vl  +  -  i  U*  Vy  +  2Uy 


V  = 

V  -  r, 


Disturbance  Continuity: 


£c£  =  p  Tdx+T  px 

£c*  =  P  T  %  +  T  py 

£C(i>  =  i  /3  pT  m) 

£cg  =  —  P  U  ft;  —  P  V  dy  —  p  Ux  —  p  Vy  —  i  5  p  W  —  U  pX  —  V  py  +  ^  {u  TX  +  V  Ty 

=  ud^+V  dy  +  \  3  VO  +  Un+Vy-^iu  Tx  +  V  Ty'j 

nci  =  -»  P 

Rep  =  1 
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